Dynamic Block Matching to assess the longitudinal component of the dense motion field of the carotid artery wall in B-mode ultrasound sequences - Association with coronary arte

Purpose: The motion of the common carotid artery tissue layers along the vessel axis during the cardiac cycle, observed in ultrasound imaging, is associated with the presence of established cardiovascular risk factors. However, the vast majority of the methods are based on the tracking of a single point, thus failing to capture the overall motion of the entire arterial wall. The aim of this work is to introduce a motion tracking framework able to simultaneously extract the trajectory of a large collection of points spanning the entire exploitable width of the image. Method: The longitudinal motion, which is the main focus of the present work, is determined in two steps. First, a series of independent block matching operations are carried out for all the tracked points. Then, an original dynamic-programming approach is exploited to regularize the collection of similarity maps and estimate the globally optimal motion over the entire vessel wall. Sixty-two atherosclerotic participants at high cardiovascular risk were involved in this study. Results: A dense displacement field, describing the longitudinal motion of the carotid far wall over time, was extracted. For each cine-loop, the method was evaluated against manual reference tracings performed on three local points, with an average absolute error of 150+/-163 um. A strong correlation was found between motion inhomogeneity and the presence of coronary artery disease (beta-coefficient=0.484, p=0.003). Conclusions: To the best of our knowledge, this is the first time that a method is specifically proposed to assess the dense motion field of the carotid far wall. This approach has potential to evaluate the (in)homogeneity of the wall dynamics. The proposed method has promising performances to improve the analysis of arterial longitudinal motion and the understanding of the underlying patho-physiological parameters.


Flow Network Tracking for Spatiotemporal and Periodic Point Matching: Applied to Cardiac Motion Analysis

The accurate quantification of left ventricular (LV) deformation/strain ...

Estimation of Cardiac Valve Annuli Motion with Deep Learning

Valve annuli motion and morphology, measured from non-invasive imaging, ...

Unsupervised Landmark Detection Based Spatiotemporal Motion Estimation for 4D Dynamic Medical Images

Motion estimation is a fundamental step in dynamic medical image process...

Left Ventricular Wall Motion Estimation by Active Polynomials for Acute Myocardial Infarction Detection

Echocardiogram (echo) is the earliest and the primary tool for identifyi...

Fast Vessel Segmentation and Tracking in Ultra High-Frequency Ultrasound Images

Ultra High Frequency Ultrasound (UHFUS) enables the visualization of hig...

Cardiac Motion Analysis by Temporal Flow Graphs

Cardiac motion analysis from B-mode ultrasound sequence is a key task in...

Early Myocardial Infarction Detection over Multi-view Echocardiography

Myocardial infarction (MI) is the leading cause of mortality in the worl...

I Introduction

Cardiovascular risk evaluation is a major public health issue as well as a tremendous scientific challenge. In the last recent years, the characterization of the elastic deformation of the tissue layers of the common carotid artery along the direction parallel to the vessel axis during the cardiac cycle in ultrasound sequences (also called cine-loops) has gathered a growing attention. This patho-physiological phenomenon, hereafter dubbed as LOKI for longitudinal kinetics, corresponds to the shear between the intima-media complex and the tunica adventitia (Figure 1persson2003new . This motion was shown to be cyclic and reproducible over several months ahlgren2012different .

Although the causes—as well as the implications—of LOKI remain to be fully determined, this dynamic parameter has been demonstrated by several groups to be associated with the presence of established cardiovascular risk factors. LOKI was shown to be associated with carotid atherosclerotic plaque burden in man and mouse svedlund2011alongitudinal , and to be a predictor for 1-year cardiovascular outcome in patients with coronary artery disease svedlund2011carotid . The motion amplitude was found to be significantly reduced in elderly diabetic patients compared to young healthy subjects zahnd2011measurement . It was also confirmed that LOKI was significantly reduced in patients with periodontal disease compared to age-matched healthy volunteers, independently of other established cardiovascular risk factors zahnd2012longitudinal . Similar results were found in a study involving individuals with spinal cord injury and able-bodied subjects, where the longitudinal retrograde motion and retrograde shear strain were reduced in the subclinical cardiovascular risk population tat2017reduced . Administration of catecholamines and in pigs showed a significant positive correlation with LOKI, thus potentially constituting a link between mental stress and cardiovascular disease ahlgren2012longitudinal . Building upon these findings, the independence between LOKI and wall shear stress was demonstrated, suggesting that LOKI is most likely caused by the pulsatile displacement of the heart, rather than by blood flow ahlgren2015profound . A similar hypothesis was formulated in a study that evaluated the association of LOKI with both left ventricular cardiac motion and local blood velocity au2016carotid . A strong association between LOKI and the severity of carotid stenosis was demonstrated, potentially constituting a tool to screen arterial sites known to be predisposed to atherosclerotic plaque formation soleimani2015novel . In addition to motion amplitude, two new indices of arterial stiffness derived from LOKI waveforms were explored: the total distance of the intima trajectory during one heart cycle, as well as the degree of a polynomial function fitting LOKI, were found to be associated with arterial stiffness measurements yli2016new . It was further demonstrated that LOKI is inversely associated with arterial stiffening taivainen2017interrelationships , and that LOKI provides complementary information regarding arteriosclerosis and risk factors compared to traditional markers taivainen2018influence . These findings therefore strongly suggest that robust and accurate LOKI quantification may improve cardiovascular risk evaluation.

Tracking the motion of the common carotid artery in vivo in B-mode ultrasound image cine-loops is hindered by several factors of different nature. First, the resolution cell of the scanner is generally coarser in the direction perpendicular to the beam () than in the direction of the beam propagation (). This is reflected by the anisotropic shape of the point spread function, whose width (determined by the probe geometry and the depth of the focal zone) is typically larger than its height (determined by the ultrasound signal wave-length) meunier1995ultrasonic . Second, the profile of the image in the longitudinal direction is rather homogeneous and does not present clear landmarks. As opposed to the radial direction, where image intensity undergoes profound variations along the depth of the image, only a small gray level variation is present along the longitudinal profile, as visible in Figure 1. This phenomenon is known as the aperture problem, and in this case, is caused by the geometry of the anatomical structure of the vessel, consisting of tissue layers aligned along the longitudinal axis. As a result, motion along the direction is more difficult to perceive. Third, images are often corrupted by several degrading phenomena inherent to the ultrasound modality itself, such as out-of-plane movements, speckle decorrelation, acoustic shadowing, or movement artifacts, which can lead to the alteration of the tracked speckle pattern over time. Another aspect to consider, especially in the context of clinical routine, is the high variability of image quality caused by the specific vessel geometry and tissue echogenicity of each examined subject.

Figure 1: (a) Longitudinal view of a common carotid artery in B-mode ultrasound. (b) Detail of the far wall. The cyclic motion of the tissue layers along the axis of the vessel (i.e., LOKI) is indicated by the double arrow.

A number of motion tracking methods have been put forward to address the aforementioned issues and evaluate LOKI. A block matching approach—also called speckle tracking—could successfully extract the longitudinal trajectory of a selected point in a given cine-loop golemati2003carotid . The motion was assessed between pairs of consecutive frames by locating the point centered on the speckle pattern bearing the greatest similarity with the tracked reference pattern, according to the sum of absolute difference metric. This study however has two limitations: First, the rather large size of the tracked blocks (

, added to the fact that the images were not interpolated, resulted in coarse, stairway-like motion trajectories. Second, this method was applied on 5 healthy subjects and 2 atherosclerotic patients, but the tracking accuracy was not quantitatively evaluated against any reference. Another approach exploited a very small block size (

and reached a high tracking accuracy when evaluated in vitro in an agar phantom with tiny glass beads producing distinct echoes cinthio2005evaluation . Reference motion was provided using laser measurements and the mean tracking error was . Nevertheless, in vivo experiments were only conducted on one healthy subject and not evaluated quantitatively. Moreover, this framework requires optimal image quality and a careful selection of a single and very well defined scatter point, thus making it not truly applicable on images acquired in routine clinical practice. A different approach was proposed based on cross-correlation after compensation for luminance changes yli2013axial . Applied in 19 healthy subjects, the evaluation protocol did not involve characterizing the accuracy of the method per se, but rather focused on quantifying the intra-subject reproducibility of the measurements on image data acquired on two subsequent days, which resulted in a Cronbach coefficient between and

. Several advanced methods based on Kalman filtering have been specifically introduced in an attempt to improve the tracking performances. An approach predicting the motion of the target point and the appearance of the corresponding speckle pattern was put forward 

gastounioti2013carotid . Although quantitative evaluation was only performed on simulated sequences, showing a 47% accuracy increase compared to traditional block matching, the method was applied on 9 healthy subjects and 31 atherosclerotic patients, demonstrating substantial difference in the measured motion amplitude. Another approach proposed the use of a control signal designed to avoid the accumulation of successive tracking errors zahnd2013evaluation . Quantitatively evaluated in 57 healthy subjects and 25 patients at high cardiovascular risk against manually traced motion references, the average accuracy () was  . A statistically significant difference was also found in the motion amplitude between controls and patients (  vs  , ). The combination of update and prediction has also been exploited within a state-space approach based on an elasticity-model and ruled by a filter gao2017robust . A quantitative evaluation of the tracking accuracy on 37 healthy subjects and 103 patients at high cardiovascular risk was performed against manually tracked points, the Pearson’s correlation coefficient was and the width of the confidence interval was .

All the methods mentioned above share a similar property: They are based on the local estimation of the temporal trajectory of a single point. This approach requires a careful selection of the best point to be processed within the image, typically, a point located in a region with optimal image quality, and centered on a well distinguishable speckle pattern. Such local analysis is limited by three principal drawbacks: First, interrogation of a unique region fails to capture the non-uniform behaviour of the dynamics of the wall, since different regions may undergo deformations of different magnitudes zahnd2015progressive . Second, tracking a single target point is likely to yield erroneous results in noisy image regions. Third, point selection—and hence results exploitation—is subject to inter-analyst variability. Therefore, characterizing the motion of the entire wall instead of a single point, which is the aim of the present method, would be a more advantageous approach to address these issues.

A few methods only have been proposed to assess the global deformation of the arterial tissues across an elongated region of interest (ROI) covering (part of) the full width of the image. The traditional block matching framework has been extended by independently tracking a small set of points and averaging the resulting trajectories zahnd2012longitudinal ; tat2017reduced

. These simple models, however, do not exploit any assumption on the physical interrelations between the motion of the different points. Individual tracking failures, although their magnitude is attenuated by the averaging operation, are therefore likely to have an impact on the resulting overall trajectory. The Velocity Vector Imaging commercial software (VVI, Research Arena 2; TomTec imaging systems GmbH, Unterschleissheim, Germany) has been applied for LOKI quantification 

svedlund2011longitudinal . However, that tool was originally designed to measure the heart dynamics in echocardiography, and the evaluation of its performances on the carotid artery showed poor results zahnd2013evaluation . Another approach based on feature detectors and descriptors to automatically identify and track several keypoints was recently introduced scaramuzzino2017longitudinal . However, three limitations can be observed: First, due to substantial computation time, this method relies on a limited amount of points (typically, 10) and therefore only captures a rather sparse distribution of the wall dynamics. Second, the coordinates of the reference points to be tracked must be reselected at each time step, with no guarantee that the same set of points will be continuously analyzed over time. Therefore the motion trajectory of a specific sub-region over time can not be captured. Third, for each pair of consecutive frames, a single displacement value is computed with the median of the motion resulting from all the tracked points. This yields a unique trajectory that describes the global displacement of the entire wall, but does not permit to study the specific displacement of different sub-regions of the wall individually.

The aim of this study is to propose an approach capable of assessing the dense motion field of the tissue layers, namely, tracking one point in every column of the image. The rationale is an attempt to minimize individual tracking errors by exploiting the collaboration between the ensemble of points. This is achieved by simultaneously extracting, in each frame of the cine-loop, the displacement of all the points by means of combinatorial analysis. More precisely, the methodological contribution of the present study is the following. The framework is based on two main steps. First, for each point (one per column), a block matching operation is conducted, exploiting three specific features (radial motion guidance with anatomical contour segmentation; use of the initial speckle pattern as a control signal to generate the tracked reference block; storage of the similarity matching describing the ensemble of potential displacements). Second, the matching score of all points are gathered in an artificial space to generate a single cost function. Then, a dynamic programming scheme is run to determine the minimal-cost path that describes the globally optimal displacement of all the points simultaneously. The framework is applied to in vivo data, validated against manual reference annotations, and assessed in the context of a pilot clinical study.

Ii Materials and Methods

The method introduced here is devised to simultaneously extract the motion of all the points located on a line. The main steps of the method—hereafter referred to as Dynamic Block Matching (DBM)—are illustrated in Figure 2. The framework is composed of two principal operations: First, a block matching search is conducted independently for all the assessed points, then the resulting similarity maps are exploited in the context of a dynamic programming scheme to determine the global optimal displacement map of the ensemble of points. The remainder of this section is organized as follow: First the image acquisition procedure is defined, then the DBM framework is presented in detail, next the evaluation protocol is described, finally a quantitative analysis between LOKI-derived parameters and cardiovascular risk factors is proposed.

Figure 2: Illustration of the main steps of the framework.

ii.1 Data acquisition and study sample

Ultrasound image acquisition was performed at the National Cerebral and Cardiovascular Center, Osaka, Japan, by a skillful medical doctor (KS), as part of the routine procedure, using a clinical scanner (Toshiba Applio, Tokyo, Japan) equipped with a linear probe (PLT-704SBT, ). Pixel size was and frame rate was 32 fps. Images were recorded in the longitudinal plane during several cardiac cycles (average number of frames in each cine-loop: , corresponding to ). The probe was oriented such that the left and right sides of the image corresponded to the heart and head direction, respectively.

Sixty-two participants suffering from atherosclerosis (age: y.o., 50 male) were enrolled in this study. The participants characteristics is provided in Table 1. Opt-out consent was obtained from all participants. This observational study was conducted as a retrospective analysis and fulfilled the requirements of our institutional review board and the ethics committee.

Parameter Unit Total Missing
Coronary artery disease 21 (34%)
Hypertension 53 (85%)
Systolic blood pressure mm Hg 14
Diastolic blood pressure mm Hg 14
Pulse pressure mm Hg 14
Dyslipidemia 55 (89%)
Peak systolic value m/s 13
Other heart diseases 27 (44%)
Diabetes mellitus 24 (39%)
Gender male, 50 (81%)
Age years 1
Previous stroke 15 (24%)
Intima-media thickness mm 8
Peripheral artery disease 8 (16%)
Estimated glomerular filtration rate mL/min/1.73m2 16
Total cholesterol mmol/L 24
Glycosylated hemoglobin % 29
Smoking: 4
- Current 7 (11%)
- Ex 39 (63%)
- Never 12 (19%)
Body mass index kg/m2 8
Other vascular diseases 6 (10%)
Table 1: Characteristics of the 62 involved participants.

ii.2 Motion tracking

ii.2.1 Initialization

First, a ROI is determined by manually indicating the left and right borders of the full exploitable width of the far wall in order to clip out noisy regions (Figure 3a). All the subsequent steps of the framework are fully-automatic.

In each frame of the image cine-loop, the contours of the lumen-intima (LI) and media-adventitia (MA) interfaces are segmented within the ROI (Figure 3a), using a previously validated method zahnd2017fully . Briefly, the contour segmentation framework is based on a front propagation approach specifically devised to simultaneously (as opposed to one after the other) extract the two anatomical interfaces of the intima-media complex. This is achieved by building a 3D space (LI depth MA depth width), where the minimum-cost path is a medial axis that describes, in each column of the image, i) the center of the intima-media complex, and ii) its thickness. The LI and MA contours can immediately be deduced from the medial axis.

In preparation for a block matching scheme, the contour information is used to populate a mesh of points , with and , in the first frame of the cine-loop. One point is placed in each column of the ROI, in the center of the intima-media complex (Figure 3b,e).

The image is then spatially interpolated: It has been shown that LOKI magnitude was approximately  cinthio2005evaluation ; golemati2003carotid ; zahnd2011measurement . Given an isotropic pixel size of , the total amplitude of the tissue motion is around 14 pixels over the cardiac cycle, which generally corresponds to 30 images. Although tissue dynamics greatly depends on the cardiac phase, the motion amplitude between two consecutive frames generally corresponds to a few pixels and may often be smaller than one pixel. Therefore, to capture the motion in finer details, a spline interpolation with a factor of is performed along the direction for sub-pixel resolution, as previously proposed by other studies cinthio2005evaluation ; zahnd2011measurement ; gastounioti2013carotid ; tat2017reduced ; gao2017robust . This operation is applied to the image, as well as the points and contour coordinates.

Figure 3: Schematic representation of the main steps of the framework. (a) Carotid ultrasound image. The lumen-intima (LI) and media-adventitia (MA) contours are segmented (yellow solid lines) within the region of interest of width (white dashed lines). The white-filled rectangle represents a block (size ), centered around a point to be tracked. (b) Detail of the intima-media complex (gray region), with the tracked block (white rectangle) centered around (red point, size pixel). (c) The search window (yellow rectangle, size  pixels in the interpolated image), yielding a 2D SSD map (same yellow rectangle), is automatically positioned in the center of the intima-media complex in the next frame. The radial displacement of is indicated by the vertical arrow. (d) Similarity vector (size ), generated from the column-wise minimum of the 2D SSD map in panel c. (e) Detail of a sub-sample of the tracked points (red points). (f) Cost function , resulting from the aggregation of the individual similarity vectors (panel d). Here (in f-h), the axis corresponds to the left-to-right indices of the points. (g) Cost functions  (front) and  (back), enforcing a unidirectional motion field (towards positive or negative values of , respectively). (h) Result of the front propagation in ruled by Equation 2, showing the optimal path (blue points). A more detailed representation of this panel is presented in Figure 4. (i) Motion field resulting from the optimal path (arrows between the red and blue points), yielding the points coordinate in the frame (in this example, ).

ii.2.2 Extraction of the raw motion

This operation consists in estimating the motion field , namely the relative -wise displacement of each point in each frame of the cine-loop. This process being carried out independently for each point, the description of the operation will be given for a single point at a given time step. The previously described interpolation scheme is applied to the frames and . A block matching operation is then carried out between these two consecutive frames: a rectangular block of size is centered around  in the  frame, and a search window of size is explored in the  frame. Matching similarity is evaluated with the sum of squared differences (SSD) metric.

Let us briefly recall three specificities that are most generally followed by traditional block matching approaches: i) the next position of the tracked point is blindly determined by the coordinates of the global minimum in the SSD map, ii) the search window in the  frame is centered around , and iii) the reference block corresponds to the speckle pattern centered around the currently tracked point . In comparison with traditional implementations, these three specific points are addressed differently by the present framework, as described hereafter.

First, the search is guided along the direction using a previously validated contour segmentation method zahnd2017fully . The rationale is to perform a first estimate of the wall radial displacement to guide motion tracking and systematically ensure that the tracked points remain within the intima-media complex. Guidance is realized by automatically setting the coordinate of the center of the search window in the center of the intima-media complex in the frame (Figure 3c). The segmentation method has been reported to delineate the LI and MA interfaces with a mean absolute error () of   and  , respectively zahnd2017fully . Therefore, block matching is conducted along the direction with a reduced maximum displacement , in order to allow a fine search around the initial position provided by contour segmentation. In contrast, since no a priori guidance is provided along the direction, a wider exploration range is permitted, with . This value has been selected based on the maximum expected displacement between two time steps and validated in a previous study zahnd2013evaluation . Given a pixel size of and an interpolation factor of in the longitudinal and radial directions, respectively, the size of the search window in pixels is equal to .

Second, the reference block used within the block matching operation is systematically updated to follow the gray-level variations of the moving images over time, while still preserving the appearance of the initial tracked pattern. The motivation of this implementation is to avoid a potential and irreversible divergence of the trajectory that may be caused by the accumulation of tracking errors due to artifacts and speckle decorrelation. This is accomplished by using a control signal , corresponding to the initial image pattern in the first frame of the cine-loop, to systematically keep the memory of the initially tracked point. To cope with small gray-level variations, the reference block is thus generated by a weighted sum of the current speckle pattern and the control signal , according to the following relation:


Here, is a weighting factor that determines the relative influence of the initial and current speckle patterns in generating the reference block.

Third, the optimal position of the ensemble of all blocks is determined simultaneously by means of a dynamic programming scheme. Please note that the aim of this paragraph is to describe the operations undertaken to prepare the dynamic programming scheme, whereas the dynamic programming algorithm per se is described thereafter in Section II.2.3. For a given tracked point , as part of the block matching operation, a search window is explored, resulting in a 2D SSD map  of size (Figure 3c). From this map, a 1D similarity vector of length is generated by selecting the smallest value in each column of the initial 2D SSD map (Figure 3d). Therefore, at this stage, no decision is taken regarding determining the displacement a single point: instead of selecting the candidate point with the lowest value in the 2D SSD map, the matching potential along the direction is stored for further use. The rationale is to simultaneously determine the optimal -wise displacement of the ensemble of all blocks by using the collection of similarity vectors in a dynamic programming scheme, as described in the following paragraph. Let us also remember that the displacement along the axis has already been determined by the contour segmentation.

ii.2.3 Computation of the motion field

In this step, at a given time , the ensemble of 1D similarity vectors , resulting from the raw motion estimation of the points , are considered collectively. A combinatorial analysis is performed to determine, among all possible solutions describing the displacement of the ensemble of all points, the optimal motion field given the following rules: i) data similarity: displacement is guided by the SSD matching criterion represented by the 1D similarity vectors , ii) non-crossing trajectories: the initial ordering of the points along the axis does not change, iii) motion smoothness: increase or decrease of the -wise distance between two neighboring points is penalized, and iv) motion uniformity: between two successive frames, all moving points must either follow a displacement towards the right or left side of the image. The implementation enforcing these rules is described below.

A dynamic programming algorithm based on front propagation is proposed. A cost function  is constructed in the artificial space , where each 1D similarity vector is oriented along the direction and centered on (Figure 3f). Here, the axis of the coordinate system corresponds to the left-to-right ranking of the points: two neighbors points and yield . The values of  outside of the region covered by the maps correspond to non-reachable positions, and are therefore set to infinity.

The condition of motion uniformity is addressed by using to generate two cost functions  (Figure 3g). In each row of (respectively, ), the value of the points whose -coordinates are strictly smaller (respectively, larger) to is set to infinity, thus enforcing a global displacement to the right (respectively, left).

A front propagation is then run to build two cumulative cost functions and . For , is initialized with the corresponding values of . Then, for , is iteratively generated according to Equation 2.


To provide a more intuitive grasp of this relation, the cumulated cost of the current node (first line, left part) is determined as the minimal value, across a reachable neighborhood , of the addition of i) the cumulated cost of the previous node (first line, right part), and ii) the sum of the cost of the current and previous nodes (second line) multiplied by a factor to penalize non-realistic distances between the two nodes (third line). Here, prevents crossing point trajectories, namely to respect the condition given two points and such that . The elasticity of the mesh is controlled by the smoothing coefficients , , and , as well as by the parameter , which represents the expected -wise distance between two adjacent points, namely one pixel in the original image corresponding to ten pixels in the interpolated image. The distance penalty is null when two neighbors points are separated by exactly ten pixels (one pixel in the non-interpolated image), and gradually becomes greater when the distance between these two points increases (stretching) or decreases (compression). A schematic representation of the front propagation is provided in Figure 4.

Figure 4: Front-propagation scheme used to generate the cumulative cost functions from the similarity vectors . For the sake of clarity, this example depicts a non-interpolated image (). Here, the direction of propagation is top-to-bottom. The gray-filled nodes depict the regions that already have been computed, and the white-filled nodes represent the regions yet to be computed. The node that is currently processed is indicated by an asterisk, and its potential neighbors are circled in black. In this example, (total number of tracked points), (reachable neighborhood), (width of the block matching window). The black lines connecting the nodes represent the successive back-tracking steps from any given node to the top of the map. Every path connecting the dots corresponds to the potential positioning of a set of points arranged along strictly increasing values. The final solution is determined as the path with the minimal cumulative cost.

Finally, the optimal paths and are extracted in both and via back-tracking zahnd2017fully from the point having the minimal cumulated cost (Figure 3h). The solution with the minimal cost among and is used to determine the motion field such that (Figure 3i).

ii.3 Method evaluation

To evaluate the performance of the method on a fair basis, a training set and a testing set were generated. The training set was composed of 20 randomly selected participants, and was used during the development phase of the method to empirically determine the optimal parameter settings. The testing set was composed of the 42 remaining participants, and was used solely during the evaluation phase of the method, using the previously determined optimal parameter settings, to determine the performance of the proposed approach on independent test data.

In order to quantify the accuracy of the method despite the lack of ground truth, manual annotations were performed to constitute a reference. For each cine-loop, a set of three reference points , were first manually placed by an experienced analyst (GZ) in the left, center, and right part of the ROI in the first frame. These anchor points corresponded to well contrasted and distinguishable speckle patterns that remained visible through the entire duration of the cine-loop. Then, the 2D trajectories of these points were manually traced during the entire cine-loop by tracking their position in all frames of the sequence. Intra-analyst variability was also assessed for all tracked points in the testing set. Manual annotations were conducted using a graphical interface that was specifically developed for this purpose, and a Wacom Intuos Pro pen tablet for improved precision.

The DBM method was then applied to the entire ROI. The three columns specifically describing the motion of the reference points were selected for evaluation. Finally, a state-of-the-art tracking method based on Kalman Block Matching (KBM) zahnd2013evaluation was used for comparison purpose and applied to track the same set of reference points . The tracking accuracy, for both the DBM and KBM methods, was evaluated by means of a point-to-point distance comparison between the resulting estimated trajectory and the reference manual tracings.

ii.4 Statistical analysis

The peak-to-peak amplitude of the motion trajectory (Figure 5c), an established LOKI-derived parameter that was demonstrated to be associated with cardiovascular risk factors zahnd2011measurement ; zahnd2012longitudinal ; tat2017reduced ; ahlgren2012longitudinal ; ahlgren2015profound ; au2016carotid , was extracted from the resulting motion field (Figure 5b,c). For each participant, was automatically determined as the average between all the local peak-to-peak amplitude values, in all the available cardiac cycles and across the entire width of the ROI.

The intra-participant variability index was introduced to evaluate the reliability of measurements based on a single point, namely to quantify LOKI reproducibility with respect to the selection of the horizontal coordinate of the initial measurement point. This parameter was assessed by means of the following procedure. For each participant, the local amplitude was automatically measured in each column of the ROI. The intra-participant variability was then established by calculating the mean absolute difference between all pairs , with , , , and the total width of the ROI.

The inhomogeneity index was introduced to quantify the consistency (or lack thereof) of the motion field across the horizontal direction of the ROI (Figure 5b). The parameter is determined as follow: A linear normalization is first applied to the motion field to scale the values between 0 and 1. The aim of this step is to generate an index that is independent of the mere amplitude of the motion. Then

is obtained by calculating the average standard deviation of the normalized motion field along the horizontal direction. The purpose of this step is to detect, in any given frame of the cine-loop, the presence of spatial regions where LOKI magnitude is greater than in other regions. Such behaviour, caused by an inhomogeneous motion field, is thus reflected by a large


Linear regression analysis was performed to evaluate the association of participants characteristics (Table 


, predictor variables) with either


(response variables). The analysis was performed separately for each response variable. A univariate model was first generated. A multivariable model was then generated, with a selection scheme of predictor variables based on three criteria:

i) association to the response variable with , as derived from the univariate model, ii) no pairwise correlation between the predictor variables (in case of two candidates variables being correlated, the one with the strongest association to the response variable was selected), and iii) adoption of the “one in twenty” rule austin2015number , leading to the selection of a total of three predictor variables. All statistical analysis was undertaken with R r2008language . The value was considered to indicate a statistically significant difference.

Iii Results

iii.1 Parameter settings

For all the 62 involved cine-loops, the average width of the processed region, corresponding the to full exploitable width of the image, was . This corresponded to an average number of tracked points , defined by the number of columns within the ROI, of . The average distance between two neighboring points used to generate the manual references was .

The optimal value for each parameter was empirically determined using the training set () by successively running the method using different parameters values and selecting the configuration that yielded the lowest tracking error. The allowed range of parameters is detailed in Table 2. The method was then applied once to the testing set (). The following parameter settings were used: size of the block , ; size of the search window, ; weight of the control signal, ; smoothing coefficients, , , .

() 0.15, 0.25, 0.30, 0.35, 0.40, 0.50, 0.55, 0.60, 0.65, 0.70
0.15, 0.25, 0.35, 0.75
0.10, 0.25, 0.50, 1.00
0.10, 0.25, 0.50
1.00, 2.00
Table 2: Allowed range of the parameter settings during the empirical training phase. Parameters with more tentative values are those having the largest impact on the performances.

iii.2 Tracking performances

The motion was successfully extracted from all the analyzed cine-loops. The quantitative evaluation of the proposed DBM framework, showing the mean absolute () difference between the trajectories resulting from the manual reference annotations and those resulting from the method, is presented in Table 3, together with the comparison against the previously introduced KBM method. Overall, the temporal trajectories resulting from the DBM method demonstrated a good similarity with the manual reference annotations, as displayed in Figure 5c. It was also observed that the DBM method performed quite robustly, as opposed to the KBM method, which in some cases failed to capture the exact tissue motion due to progressive drift or sudden jumps to another target point. This is reflected by the Figure 6

, where outliers are more numerous in the case of the KBM method.

Figure 5: Representative example result from the testing set. (a) Image depicting the region of interest of width (ROI, white dashes) and the set of three points used for evaluation (red circles). (b) Dense LOKI field (2D) resulting from the DBM method, displaying the longitudinal motion across the entire width of the ROI. The red dashes indicate the trajectory of the points . (c) Trajectories (1D) of the three evaluated points, for the DBM method (purple), the manual reference tracings (green), and the KBM method (orange). The tracking failures of the KBM method are indicated by the arrows. The peak-to-peak amplitude of the motion trajectory is represented. The blue dots on the temporal axes (b,c) represent the end-diastole.
Method Training set () Testing set ()
Table 3: Average absolute tracking error ( standard deviation) in compared against manual reference tracings for the Dynamic Block Matching (DBM) and Kalman Block Matching (KBM) methods, and intra-analyst variability (MAN).

The very nature of the DBM framework enabled for the extraction of a dense motion field . Several dynamic patterns, corresponding to different participants, are displayed in Figure 7. Although a quantitative analysis of the tracking accuracy was only performed in three local points per cine-loop, as described above, a qualitative assessment of the motion homogeneity can be obtained from these 2D fields.

Figure 6: Illustration of the absolute tracking errors, for all the 42 participants of the testing base, for the DBM (purple) and KBM (orange) methods. Percentiles are indicated by boxes (25th and 75th), inner lines (50th) and error bars (5th and 95th). The level is represented by the black horizontal line.
Figure 7: Example result of the longitudinal motion field during several cardiac cycles, evaluated in the intima-media complex of the far wall (arrow in panel a) within the region of interest of width (white dashes), for six different participants of the testing set. (a-c) The motion field is quite homogeneous. (d-f) The motion field is rather inhomogeneous, substantial variations are visible in both spatial and temporal directions. The motion amplitude in is indicated by the colorbars.

The DBM framework was implemented in MATLAB (MATLAB R2016b, The MathWorks Inc., Natick, MA, USA, 2016) and accelerated using Mex/C++, on a desktop computer with a 2.90 GHz processor and 32 Gb RAM. The required computational time to process a typical cine-loop (151 frames, width , corresponding to tracked points) was . More specifically, it was per frame: for the layers segmentation, for the block matching operation, and for the dynamic programming operation.

iii.3 Intra-participant variability

The intra-participant variability of LOKI with respect to the selection of the -coordinate of the initial measurement point was assessed by means of the index, as detailed in Section II.4. Among all participants, the mean absolute () value of was . With an average peak-to-peak amplitude of , the average relative intra-participant variability was of the total motion amplitude.

iii.4 Motion inhomogeneity

The association between either the motion inhomogeneity index or the peak-to-peak amplitude and participants characteristics is provided in Table 4. The strongest correlates to were coronary artery disease (, , Figure 8) and systolic blood pressure (, ). Interestingly, the association with coronary artery disease remained nearly-significant (, ) in a multivariable model after adjusting for systolic blood pressure and peak systolic value, the second and third strongest independent correlates. The peak-to-peak amplitude was strongly correlated with peripheral artery disease (, ) and hypertension (, ). These results are in accordance with previous findings showing that is associated with cardiovascular risk factors zahnd2012longitudinal . A significant association was also observed between and (, ). A qualitative analysis of the motion field in all participants did not reveal an association between motion homogeneity and putative local tissue stiffness (e.g. healthy sections or atherosclerotic plaques, Figure 7).

Inhomogeneity index Peak-to-peak amplitude
Univariate Multivariable Univariate Multivariable
Coronary artery disease .586 .003 (**) .484 .061 (.) .000 .995
Hypertension .577 .033 (*) .120 .040 (*) .094 .151
Systolic blood pressure .015 .053 (.) .008 .357 -.003 .151
Diastolic blood pressure .071 .552 -.002 .532
Pulse pressure .015 .080 (.) -.002 .240
Dyslipidemia .496 .102 -.006 .923
Peak systolic value -.003 .111 -.002 .253 .000 .198
Other heart diseases .293 .131 .054 .196
Diabetes mellitus .206 .300 -.018 .683
Gender -.235 .336 .042 .430
Age .014 .372 .000 .898
Previous stroke .162 .473 -.018 .719
Intima-media thickness -.080 .520 -.002 .936
Peripheral artery disease .154 .595 .165 .006 (**) .149 .012 (*)
Estimated glomerular filtration rate -.003 .602 -.001 .447
Total cholesterol -.002 .699 .001 .153
Glycosylated hemoglobin -.050 .725 .006 .820
- Current -.130 .773 .017 .857
- Ex .107 .770 -.072 .354
- Never -.134 .747 .006 .949
Body mass index .006 .863 -.012 .079 (.) -.012 .069 (.)
Other vascular diseases .009 .976 .047 .507
Table 4: Association of cardiovascular risk factors with LOKI-derived parameters in all 62 participants. The symbols describing the value are 0 (***) 0.001 (**) 0.01 (*) 0.05 (.) 0.1 ( ) 1.
Figure 8: Univariate association between LOKI inhomogeneity index and the presence of coronary artery disease (CAD). Percentiles are indicated by boxes (25th and 75th), inner lines (50th), and error bars (5th and 95th). The result of the Mann-Whitney U-test is indicated by the value.

Iv Discussion

The introduced method is devised to simultaneously analyze the motion of a large number of tracked points arranged in a linear fashion. It was applied in vivo in carotid B-mode ultrasound cine-loops to extract LOKI, the tissue motion along the vessel axis during the cardiac cycle. The motion is extracted along the entire exploitable length of the imaged artery (as opposed to traditional approaches that only focus on a single point). Such dense motion field analysis was a missing link in LOKI-based studies.

Let us start by briefly recalling the three hypotheses that were investigated in the present study. First, tracking of a single point is more prone to generate erroneous trajectories compared to exploiting the collaboration of multiple points. The increased robustness of DBM was demonstrated in Section III.2 against a state-of-the-art method based on single-point tracking, with an average absolute error of compared to . Second, the selection of a single point is an analyst-dependent procedure and may hinder the reliability of results exploitation. The impact caused by single-point selection when measuring the peak-to-peak motion amplitude was assessed in Section III.3, and showed that the intra-participant variability was . Although relatively reduced compared to the total magnitude of the motion, this additional source of error can be avoided with a global analysis. Third, motion trajectories resulting from a single point do not allow the analysis of the displacement field (in)homogeneity across the entire length of the artery, which may reflect cardiovascular health. The results presented in Section III.4 establish that the LOKI-derived inhomogeneity index is strongly associated with coronary heart disease (, , Figure 8).

The collaboration of all the points allowed for improved robustness. The spatial consistency of the motion was implicitly exploited within the front propagation step (Eq. 2), since the global optimum was preferred over a series of independent best matches, thus reducing the risk of individual tracking errors. To phrase it differently, the set of conditions ruling the behaviour between neighboring points could successfully maintain the cohesion of the mesh and prevent sudden jumps or cumulative drift errors. On the other hand, the mesh was flexible enough for the points to follow the motion of the tissues over time. This phenomenon is illustrated in Figure 5c. At this point, it is useful to differentiate accuracy and robustness. While accuracy describes the ability of a method to yield results close to the ground truth, robustness represents the capacity of a method to avoid failure and provide an exploitable result. In the present situation, it is unlikely that the DBM method could have a positive impact on the tracking accuracy per se, however the robustness was improved. Moreover, it should be noted that the combinatorial analysis, solved by means of dynamic programming, also enabled a simultaneous motion computation for a large collection of points in a rather limited time due to the intrinsic reduction of the solution space.

The motion along the direction was guided by the contour segmentation and a reduced search window, whereas the motion along the direction, substantially more challenging to assess, was obtained by means of a dynamic programming approach. The user interaction was minimal, only the left and right borders of the ROI to be processed had to be selected to discard non-exploitable regions of the image. In addition to providing motion information for the entire ROI, this procedure (namely, region selection rather than point selection) reduces user variability compared to individual point selection and facilitates full-automation.

A fair evaluation was conducted by using a training set () to define the optimal parameter settings and an independent testing set () to quantify the accuracy of the method. Results demonstrated that the DBM framework was able to perform in the testing set almost as well as in the training set, showing the ability of the method to analyze previously unseen data. It is also noteworthy that the method was evaluated in challenging conditions, as image quality was overall quite poor, mostly for two reasons. First, images acquired with the ultrasound scanner were rather coarse (pixel size of ). Second, the involved participants were all elderly patients suffering from atherosclerosis, and therefore data was more challenging to process than images from young healthy volunteers due to a less favorable echogenicity.

Although the motion inhomogeneity index is amplitude-independent since it is derived from the normalized motion field, this index was correlated to the peak-to-peak amplitude (, ). Nevertheless, confronting those parameters to the participants characteristics indicated that and are strongly correlated to different predictor variables (Table 4). These preliminary findings suggest that do not simply replicate the information already captured by but instead represents a complementary measure of cardiovascular risk. A comprehensive analysis of the association of motion-derived parameters with risk factors was beyond the scope of the present study. Future work involving a greater number of participants will be conducted to confirm the clinical usefulness of analyzing the spatial inhomogeneity of the motion along the wall to characterize at-risk individuals.

V Conclusion

The main contribution of this article is the introduction of a novel motion estimation framework, called DBM for Dynamic Block Matching, specifically designed to extract the temporal trajectory of a large collection of points along a direction of interest. The method was applied in ultrasound cine-loops of the common carotid artery to analyze the motion of the intima-media complex in the direction parallel to the blood flow (i.e., longitudinal kinetics, LOKI) and extract a dense motion field. Evaluated in vivo in 42 participants, the proposed method demonstrated good tracking performances compared to manually established references. The analysis of the inhomogeneity of the tissue motion across the entire length of the artery showed a strong association with the presence of coronary artery disease. These findings suggest that the proposed method may provide complementary information for cardiovascular risk assessment.

Conflict of interest

No benefits in any form have been or will be received from a commercial party related directly or indirectly to the subject of this manuscript.

This work was partly supported by the JSPS #PE16208 funding as well as the MEXT KAKANHI #26108004 funding.


  • (1) Persson M, Rydén-Ahlgren Å, Jansson T, Eriksson A, Persson HW, Lindström K. A new non-invasive ultrasonic method for simultaneous measurements of longitudinal and radial arterial wall movements: first in vivo trial. Clin Physiol Funct Imaging. 2003;23(5):247–251.
  • (2) Rydén-Ahlgren Å, Cinthio M, Persson HW, Lindström K. Different patterns of longitudinal displacement of the common carotid artery wall in healthy humans are stable over a four-month period. Ultrasound Med Biol. 2012;38(6):916–925.
  • (3) Svedlund S, Gan L. Longitudinal common carotid artery wall motion is associated with plaque burden in man and mouse. Atherosclerosis. 2011;217(1):120–124.
  • (4) Svedlund S, Eklund C, Robertsson P, Lomsky M, Gan LM. Carotid artery longitudinal displacement predicts 1-year cardiovascular outcome in patients with suspected coronary artery disease. Arterioscler Thromb Vasc Biol. 2011;31(7):1668–1674.
  • (5) Zahnd G, Boussel L, Marion A, et al. Measurement of two-dimensional movement parameters of the carotid artery wall for early detection of arteriosclerosis: a preliminary clinical study. Ultrasound Med Biol. 2011;37(9):1421–1429.
  • (6) Zahnd G, Vray D, Sérusclat A, et al. Longitudinal displacement of the carotid wall and cardiovascular risk factors: associations with aging, adiposity, blood pressure and periodontal disease independent of cross-sectional distensibility and intima-media thickness. Ultrasound Med Biol. 2012;38(10):1705–1715.
  • (7) Tat J, Au JS, Keir PJ, MacDonald MJ. Reduced common carotid artery longitudinal wall motion and intramural shear strain in individuals with elevated cardiovascular disease risk using speckle tracking. Clin Physiol Funct Imaging. 2017;37(2):106–116.
  • (8) Rydén-Ahlgren Å, Cinthio M, Steen S, et al. Longitudinal displacement and intramural shear strain of the porcine carotid artery undergo profound changes in response to catecholamines. Am J Physiol Heart Circ Physiol. 2012;302(5):H1102–H1115.
  • (9) Rydén-Ahlgren Å, Steen S, Segstedt S, et al. Profound increase in longitudinal displacements of the porcine carotid artery wall can take place independently of wall shear stress: a continuation report. Ultrasound Med Biol. 2015;41(5):1342–1353.
  • (10) Au JS, Ditor DS, MacDonald MJ, Stöhr EJ. Carotid artery longitudinal wall motion is associated with local blood velocity and left ventricular rotational, but not longitudinal, mechanics. Physiol Rep. 2016;4(14):e12872.
  • (11) Soleimani E, Mokhtari-Dizaji M, Saberi H. A novel non-invasive ultrasonic method to assess total axial stress of the common carotid artery wall in healthy and atherosclerotic men. J Biomech. 2015;48(10):1860–1867.
  • (12) Yli-Ollila H, Laitinen TP, Weckström M, Laitinen TM. New indices of arterial stiffness measured from longitudinal motion of common carotid artery in relation to reference methods, a pilot study. Clin Physiol Funct Imaging. 2016;36(5):376–388.
  • (13) Taivainen SH, Yli-Ollila H, Juonala M, et al. Interrelationships between indices of longitudinal movement of the common carotid artery wall and the conventional measures of subclinical arteriosclerosis. Clin Physiol Funct Imaging. 2017;37(3):305–313.
  • (14) Taivainen SH, Yli-Ollila H, Juonala M, et al. Influence of cardiovascular risk factors on longitudinal motion of the common carotid artery wall. Atherosclerosis. 2018;272:54–59.
  • (15) Meunier J, Bertrand M. Ultrasonic texture motion analysis: theory and simulation. IEEE Trans Med Imaging. 1995;14(2):293–300.
  • (16) Golemati S, Sassano A, Lever MJ, Bharath AA, Dhanjil S, Nicolaides AN. Carotid artery wall motion estimated from B-mode ultrasound using region tracking and block matching. Ultrasound Med Biol. 2003;29(3):387–399.
  • (17) Cinthio M, Rydén-Ahlgren Å, Jansson T, Eriksson A, Persson HW, Lindström K. Evaluation of an ultrasonic echo-tracking method for measurements of arterial wall movements in two dimensions. IEEE Trans Ultrason Ferroelectr Freq Control. 2005;52(8):1300–1311.
  • (18) Yli-Ollila H, Laitinen T, Weckström M, Laitinen TM. Axial and radial waveforms in common carotid artery: an advanced method for studying arterial elastic properties in ultrasound imaging. Ultrasound Med Biol. 2013;39(7):1168–1177.
  • (19) Gastounioti A, Golemati S, Stoitsis JS, Nikita KS. Carotid artery wall motion analysis from B-mode ultrasound using adaptive block matching: in silico evaluation and in vivo application. Phys Med Biol. 2013;58(24):8647.
  • (20) Zahnd G, Orkisz M, Sérusclat A, Moulin P, Vray D. Evaluation of a Kalman-based block matching method to assess the bi-dimensional motion of the carotid artery wall in B-mode ultrasound sequences. Med Image Anal. 2013;17(5):573 - 585.
  • (21) Gao Z, Xiong H, Liu X, et al. Robust estimation of carotid artery wall motion using the elasticity-based state-space approach. Med Image Anal. 2017;37:1–21.
  • (22) Zahnd G, Balocco S, Sérusclat A, Moulin P, Orkisz M, Vray D. Progressive attenuation of the longitudinal kinetics in the common carotid artery: preliminary in vivo assessment. Ultrasound Med Biol. 2015;41(1):339–345.
  • (23) Svedlund S, Gan LM. Longitudinal wall motion of the common carotid artery can be assessed by velocity vector imaging. Clin Physiol Funct Imaging. 2011;31(1):32–38.
  • (24) Scaramuzzino S, Carallo C, Pileggi G, Gnasso A, Spadea MF. Longitudinal motion assessment of the carotid artery using speckle tracking and scale-invariant feature transform. Ann Biomed Eng. 2017;45(8):1865–1876.
  • (25) Zahnd G, Kapellas K, Hattem M, et al. A fully-automatic method to segment the carotid artery layers in ultrasound imaging – Application to quantify the compression-decompression pattern of the intima-media complex during the cardiac cycle. Ultrasound Med Biol. 2017;43(1):239–257.
  • (26) Tat J, Au JS, Keir PJ, MacDonald MJ. Reduced common carotid artery longitudinal wall motion and intramural shear strain in individuals with elevated cardiovascular disease risk using speckle tracking. Clin Physiol Funct Imaging. 2017;37(2):106-116
  • (27) Austin PC, Steyerberg EW. The number of subjects per variable required in linear regression analyses. J Clin Epidemiol. 2015;68(6):627–636.
  • (28) R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2008. ISBN 3-900051-07-0.