Digital image correlation for full-field time-resolved assessment of arterial stiffness

Abstract. Pulse wave velocity (PWV) of the arterial system is a very important parameter to evaluate cardiovascular health. Currently, however, there is no golden standard for PWV measurement. Digital image correlation (DIC) was used for full-field time-resolved assessment of displacement, velocity, acceleration, and strains of the skin in the neck directly above the common carotid artery. By assessing these parameters, propagation of the pulse wave could be tracked, leading to a new method for PWV detection based on DIC. The method was tested on five healthy subjects. As a means of validation, PWV was measured with ultrasound (US) as well. Measured PWV values were between 3.68 and 5.19  m/s as measured with DIC and between 5.14 and 6.58  m/s as measured with US, with a maximum absolute difference of 2.78  m/s between the two methods. DIC measurements of the neck region can serve as a test base for determining a robust strategy for PWV detection, they can serve as reference for three-dimensional fluid–structure interaction models, or they may even evolve into a screening method of their own. Moreover, full-field, time-resolved DIC can be adapted for other applications in biomechanics.


Introduction
In digital image correlation (DIC), a speckle pattern is physically applied to an object. Coordinates of points, labeled by the randomly applied stochastic (speckled) pattern are captured with a camera and identified with an image correlation. 1 These points are followed when the object is deformed and displacements are acquired with subpixel resolution. Full-field, threedimensional (3-D) surface results can be acquired using a stereo camera system. The optical setup for 3-D acquisition consists of two cameras and a bright light source only, making this technique easy to deploy and measurements easy to perform. Although image analysis is computationally intensive, increasing computer power has made the technique more popular 2 resulting in biomechanical applications such as strain measurements on human bone and tendon, 3 on mouse arteries, 4 on mouse bone, 5 and on finch beak. 6 In the present work, DIC is applied for the first time, to our knowledge, on a living (human) subject. DIC was used for timeresolved, full-field assessment of the movement of the skin in the neck directly above the common carotid artery (CCA). The movement of the skin on this location is directly linked to the pressure inside the CCA (Ref. 7). Detecting skin motion on this spot thus allows tracking the pressure wave inside CCA. This pressure wave propagates with a certain velocity, the pulse wave velocity (PWV). PWV of the arterial system is a measure for cardiovascular health. 8 Several methods to measure PWV exist today. Currently, PWV is being measured between the femoral artery (FA) in the groin and the CCA in the neck. For this purpose, there are contact methods based on tonometry 9 or ultrasound (US), 10 or noncontact methods based on laser Doppler vibrometry (LDV). 11 The main advantage of this approach is that it is already well integrated in the clinical circuit. However, the patient is required to be undressed and PWV detection between the neck and the groin poses specific difficulties. 12 Recently, local PWVof the CCA has attracted a serious interest for screening purposes. The CCA is often affected by arteriosclerosis and atherosclerosis. 13 Also, elevated PWV in this part of the arterial system is hypothesized to play a major role in the etiology of cardiovascular disease. 14 Moreover, the easy access of the artery allows the patient to remain dressed during measurements. Local PWV measurements in the CCA have been realized using US (Refs. 15 and 16) and LDV. 17 However, this approach is relatively new in general and it lacks validation at the moment.
In this work, DIC is proposed as a screening method on its own for local PWV detection of the CCA. Data acquired with DIC describe full-field, time-resolved deformation of the skin overlaying the CCA. This information can serve as a test base for determining a robust strategy for PWV detection as they create insight in the movements of the skin above the artery and the pressures inside the artery. Moreover, they can serve as a reference for 3-D fluid-structure interaction (FSI) models of the CCA when complete systems, including arteries and surrounding tissue and skin, are being modeled. Finally, the measurement approach can be adapted easily for other applications in biomechanics such as determination or modeling of biomaterials, study of growth dynamics, among other possible examples.

Measurement Setup
A stereo system with two cameras (Redlake HR1000 and Redlake Motion PRO, IDT, Tallahassee, Florida), with 8-bit depth, resolution of 1280 × 1024 pixels, and maximum frame rate of 500 fps, was used. The cameras were placed at a distance of approximately 1 m from the subject and the mutual distance between cameras was approximately 0.5 m. Both cameras were mounted with an identical 200 mm f/4 lens (Nikon, Tokyo, Japan) with an opening angle of approximately 30 deg. Movies were synchronized post hoc, using a short flash from a laser pen. High power light-emitting diodes (LEDs) fed with distensibility coefficient (DC) to avoid flickering illuminated the sample (see Fig. 1). Automatic system calibration was performed by using various images of a translated and rotated regular grid pattern within a bundle adjustment technique taken before, between, and after the entire experiment. 18,19 As a result, intrinsic and extrinsic camera parameters were obtained (see Table 1). DIC requires the presence of a random pattern on the specimen. This was applied by air brushing black paint (Air Stream Makeup Black, Kryolan, Berlin, Germany) on top of a white ground layer (Wet Makeup White, Kryolan, Berlin, Germany) in the left neck region at the level of the left CCA. In order for DIC to perform optimally, the air brushed speckle pattern should have a random character, and the contrast between the speckles and the backgrounds should be as high as possible. Without the preparation of the surface, i.e., on the bare skin, DIC will not render satisfying results. Five young male healthy volunteers participated in this study. Each measurement, a movie of about 1 to 2 s, was recorded with a sample frequency of 500 images∕s, in order to capture the movements generated by between 3 and 6 heartbeats. (The exact number of recorded heartbeats per subject is given in Sec. 3.) The whole measurement was repeated one-day later with a different speckle pattern on the same location. In one image, the initial nondeformed shape was reconstructed by finding the corresponding speckles in the images captured by camera 1 and camera 2 and via a triangulation method invoking the determined stereo camera parameters. 20 In the next step, a similar procedure was applied to the subsequent image of the deformed state. Finally, a direct comparison of the deformed ðx; y; zÞ and the undeformed ðx0; y0; z0Þ shape yields the 3-D deformation coordinates ðu; v; wÞ. These contain all the information needed to determine the in-plane normal and shear strain components ðε xx ; ε yy ; ε xy Þ, or the principal strains and the shear angle ðε 1 ; ε 2 ; γÞ as described in Ref. 21 The following settings were used in the displacement determination: a subset size of 16 × 16 pixels that is being compared between two images, a step size of 5 pixels determining the resolution of the correlation analysis, zeronormalized cross-correlation algorithm, 22 bicubic interpolation, 23 and affine shape functions. 24 Strains were calculated in the Green-Lagrange convention 25 and a strain window of 50 × 50 pixels and a bilinear interpolation 21 was taken. The results are transformed into the coordinate system of camera 1. All calibration and correlation steps are executed in the MatchID software package (in-house developed: http://www.matchid.org/).

Data Analysis
For further analysis, the deformation in the out-of-plane direction (z) of every frame was assessed with the dedicated algorithms implemented in MATLAB (MathWorks, Natick, Massachusetts). First z-displacement was filtered with a 3-D polynomial filter. (Order was 1 × 1 × 3 and 3-D window size was 10 × 10 × 15 in the x, y, and time dimension, respectively.) The filter first provided out-of-plane displacement, but also the velocity (first-derivative) and acceleration (second-derivative) in the time domain in one filtering step (see Fig. 2). Although the Fig. 1 The experimental setup with the relative position of the stereo system and the measurement area (left panel). x and y axes define the in-plane components, out-of-plane components are reported in the right panel as the z-axis. A stereo system with two cameras with an opening angle of approximately 30 deg was used. Approximate distance D1 between cameras and subjects was 1 m, approximate mutual distance between cameras was 0.5 m. High power LEDs fed with distensibility coefficient to avoid flickering illuminate the sample. In the neck in the region around the CCA, a ground layer of white paint is applied, covered with a black speckle pattern. The trajectory of the artery is indicated with a black arrow. The white cross in the right panel refers to Fig. 2. paint was applied around the presumed trajectory of the artery as found by palpation, and although the trajectory was indicated with markings from a marker pen, arterial trajectory was estimated in an extra step by plotting maximal displacement during one heartbeat [see Figs. 3(a) and 3(b)], hypothesizing that the distension of the skin lying on top of the tissues flanking the artery (muscle and larynx) is much smaller than the distension of the skin lying on top of the artery. Then, the time domain was resampled 30 times in MATLAB using a polyphase filter implementation, 26 increasing sample frequency to 15;000 samples∕s. For further processing, the acceleration of the skin surface in the z-direction was considered (see Fig. 1). For each pixel, the Table 1 Automatic system calibration was performed by using various images of a translated and rotated regular grid pattern within a bundle adjustment technique taken before, between, and after the entire experiment. As a result, intrinsic and extrinsic camera parameters were obtained.

Intrinsic parameters
Extrinsic parameters   acceleration of the skin surface in the z-direction was analyzed for the two specific local maxima, commonly used for the detection of PWV: the maximum of z-acceleration of the surface and the dicrotic notch 27 (see Fig. 2). To detect these maxima, first, a time window was manually delineated in one randomly chosen pixel. Then the local maximum was determined automatically in this window in every other pixel. After detection of the maxima, a visual inspection of accuracy of this approach was performed. Maps were made with the relative progress of the wave (see Fig. 4). Second, the same characteristic time points were analyzed along the trajectory of the artery as described in the previous paragraph (see Fig. 3). The distance along the trajectory of the artery was plotted against the apparition in time of the local maxima, and from the slope of this relationship, the velocity of the wave, or the PWV, was derived and its measurement error (ME) (see Fig. 5). PWV was determined for every recorded pulse wave using the two described algorithms (see Tables 2  and 3). The strains were not subjected to any additional processing and are mentioned for illustrative purposes only.

Ultrasound Measurements
In order to compare the measured PWV values with the parameters of the CCA, in each test subject arterial parameters were determined with the US 1 to 3 weeks before the DIC measurements, performed by trained clinicians of the Department Fig. 4 For each pixel, time domain was analyzed for the characteristic time points in the shape of each recorded pulse wave, in order to track the wave and determine its velocity. This was done for the entire measured domain for two different characteristic time points commonly used in PWV detection: the maximum of z-acceleration of the surface (a), and the dicrotic notch (b). Contour plots were made using the relative progress of these time points, thus displaying the progress of the wave. Again, the trajectory of the artery is indicated with a black arrow. Isolines are in milliseconds. In the upper right corner of this figures, some noise is present. Fig. 5 The trajectory was found as explained above. The trajectory is displayed in Fig. 3. The distance along the trajectory of the artery was plotted against the moment of appearance of the characteristic time points, and from the slope of this relationship, the velocity of the wave, or the PWV, was derived. (a) PWV is calculated by finding the maximum of z-acceleration and (b) PWV is calculated by finding the dicrotic notch. In this particular example, we find for PWV acceleration : PWV acceleration ¼ ð4.00 AE 0.04Þ m∕s, Cardiology of the University Hospital of Antwerp. All measurements were performed three times in a row. Test subjects were asked to remain sober 3 h prior to the US measurement. Arterial parameters were used to estimate PWV, using the Bramwell-Hill equation: where A is the arterial cross-sectional area during diastole, dA is the difference in the cross-sectional area between systole and diastole, ρ is the blood density, and dP is the pressure difference in the artery between systole and diastole. 7 The US and DIC methods were compared using a Bland-Altman (BA) analysis. For comparison, the average of PWV pooled per individual, per method-and if applicable-per algorithm was used. For the BA analysis, a range of agreement was defined as mean bias AE2 times standard deviation (SD). 28 3 Results

Digital Image Correlation
Out-of-plane displacement, velocity, and acceleration of the CCA region were calculated in 3-D, visualizing pulse wave propagation in the arterial segment. When the time domain of the signal is analyzed for the specific points in the shape of the pulse wave, such as the maximum of acceleration, or the dicrotic notch, maps can be made of the relative progress of the pulse wave (see Fig. 4). When the signal was considered along the trajectory of the artery (see Fig. 3), it can be observed that the pulse wave propagates in the segment (see Fig. 6). Additionally, PWV could be calculated from these waveforms along the arterial segment (see Fig. 5). PWV acceleration of individual heartbeats was found to be between (3.17 AE 0.08) and ð5.32 AE 0.05Þ m∕s using the maximum of acceleration (reported as value AE ME), while average PWV acceleration per subject was found to be between (3.68 AE 0.78) and ð5.09AE 0.22Þ m∕s (reported as average AE SD). PWV dicrotic notch of individual heartbeats was found to be between (3.04 AE 0.30) and ð6.23 AE 0.66Þ m∕s using the dicrotic notch (reported as value AE ME), while average PWV dicrotic notch per subject was found to be between (4.05 AE 0.39) and ð5.19 AE 0.53Þ m∕s (reported as average AE SD). PWV as measured with DIC is summarized in Fig. 7 and Tables 2 and 3. Also, strains of the skin surface were calculated (see Figs. 8 and 9).

Ultrasound
Blood pressure, intima media thickness, lumen diameter, distension, compliance coefficient, DC, and PWV arterial parameters were determined (see Table 4). PWV was found to be between 5.24 AE 0.22 and 6.58 AE 0.32 as calculated from the former parameters (see Fig. 7 and Table 4). Since PWV values are calculated from averaged parameters with SD, PWV results are displayed as average with ME.
The BA analysis indicates that the 95% limits of agreement between the PWV acceleration compared to PWV arterial parameters ranged from −3.23 to 0.77 m∕s. The two methods do not consistently provide similar measures because there is a level of disagreement that includes discrepancies of up to −1.23 m∕s. Table 2 In each of the five test subjects PWV was determined twice (measurements 1 and 2) using two different algorithms (PWV acceleration and PWV dicrotic notch ). N indicates the amount of samples (heartbeats) for each (type of) measurement. Results are reported as maximum or minimum recorded value AE ME.

Description
PWV acceleration PWV dicrotic notch Maximum N AEME (m∕s) AEME (m∕s) AEME (m∕s) AEME (m∕s) AEME (m∕s) AEME (m∕s) AEME (m∕s) AEME (m∕s) The BA analysis indicates that the 95% limits of agreement between the PWV dicrotic notch compared to PWV arterial parameters ranged from −2.85 to 0.44 m∕s. The two methods do not consistently provide similar measures because there is a level of disagreement that includes discrepancies of up to −1.20 m∕s (see Fig. 10).

Discussion and Conclusion
PWV is a clinical important marker for cardiovascular disease. 8 However, a golden standard for PWV detection is nonexistent at the moment. 29 Several approaches are available. Most methods focus on the PWV as measured between the carotid and the FA.   This approach measures PWV over the entire arterial system, while stiffness of central and peripheral arteries often differs. 30 Also, pulse wave form will change dramatically as measurement spots lie further apart, making it hard to define pulse transit time. 31 Moreover, estimating the distance between measurement sites can introduce errors. 32 Finally, the patient is required to undress during measurements, which can be an issue in some parts of the world. Aforementioned reasons sparked interest in alternative ways of PWV determination. Currently, determination of local PWV in the CCA attracts scientific attention. The CCA is an interesting target: the artery surfaces at an easily accessible place, distance between measurement sites can be readily measured, pulse wave will be minimally changed over short assessed distances, it is often affected by atherosclerosis, and augmented stiffness of the CCA is believed to cause hypertrophy of the heart due to the wave reflections at the bifurcation. 33 Several noninvasive methods focus on the PWV detection in the CCA, such as pulse wave imaging (PWI), 16 and LDV. 17,34 Temporal resolution of reported PWI studies is often rather low (order of several 100 Hz), and LDV has a limited amount of measurement points (two measurement points when two Fig. 9 Full-field strains were calculated for the skin above the CCA. Here, strains and z-displacements are displayed between measurement points as indicated in Fig. 8. The left panel displays displacement, the middle panel displays strain E1, and the right panel displays strain E2. In the plane displaying displacement and in the panel displaying E1 strain, a clear longitudinal region is highlighted along the trajectory of the artery. Table 4 In each of the five test subjects arterial parameters were determined with the US. Wall thickness (IMT), lumen diameter (D), distension, and PWV arterial parameters were determined. N indicates the amount of samples (heartbeats) for each (type of) measurement. Results are given as mean AE SD, except for PWV which is given as mean AE ME as it is calculated from several parameters with the SD. The relative difference of the averaged PWV values as acquired with DIC is also given, using US as the reference.

Subject (#)
Arterial parameters Relative difference (%) regular LDV systems are used). Moreover, aiming the LDV systems exactly on top of the artery can be difficult.
In the present work, a DIC based noncontact method is presented for the local PWV measurements in the CCA.
Temporal resolution of 500 frames∕s is achieved with the current setup, equaling resolution acquired with US. For young healthy volunteers, where the PWV is low, this frame rate appears to be sufficient to detect PWV. However, in more clinically relevant populations, PWV will attain values of 13 m∕s and more. 8 As such, a higher frame rate will be required. Frame rates up to several kilohertz are within the reach of modern high-speed camera systems and will be the subject of future research in order to attain validation in a population with higher PWV. Beat-to-beat variability defined as SD is as high as 1.56 m∕s with DIC, and intersession variability is as high as 1.35 m∕s. Since local PWV measurements are not yet established as a screening method for cardiovascular disease, the question whether this is sufficient has yet to be answered. However, it can be expected that intersession and beat-to-beat variability of the measurements will improve as frame rates increase.
Spatial resolution, on the other hand, is limited by factors such as resolution of the cameras at a specific sample frequency, position of the cameras, quality of the calibration step, characteristics of the lenses, coarseness of the speckle pattern, and correlation software. 35 An experimental or theoretical determination of the uncertainties on presented quantities was not performed during this experiment, since we were only interested in the actual shape of the waveforms. However, post hoc, we were able to achieve a resolution better than 0.015 mm.
Due to the achieved spatial resolution, it was possible to visualize the trajectory of the artery (see Fig. 3) highlighting the importance of knowing the exact trajectory of the artery for local methods for PWV detection. Subsequently, PWV was assessed along this trajectory using known methods for the pulse wave detection (Fig. 5). Finally, maps of the relative progression of the wave could be made (see Fig. 4).
However, it can be seen that the isolines in Fig. 4 are nonequidistant. This could be explained by nonideal elastic behavior and structural inhomogeneity in the trajectory of the artery or by artifacts due to undersampling and smoothing of the data, or by a combination of both. Also, the isolines appear to be slightly off axis. The latter could be due to a slight tilt of the direction in which pulsation is at its maximum, and the z-axis of the measurement setup. To tackle this problem, coordinate systems could be rotated, such that the displacements in the z-direction are maximal (e.g., using principal component analysis).
In order to validate the results, PWV values obtained with DIC measurements were compared with those obtained from the US measurements as derived with the Bramwell-Hill equation (1) (see Fig. 7). However, validation of the local PWV measurements poses specific problems. The presented method is a time-of-flight approach, typically prone to errors introduced by wave reflections and changes in waveform over distance. 12,29,33 A single point measurement method could circumvent this problem. However, this would require extra information, such as local blood velocity profiles or local blood pressure profiles. 36,37 It can be expected that the additional measurement equipment for this purpose would make the technique either impractical or invasive.
PWV velocities as obtained from the Bramwell-Hill equation (1), however, give only an approximation of the "ground-truth" PWV, assuming ideal circumstances. PWV as calculated with this equation assumes the absence of wave reflections, homogeneous mechanical parameters along the arterial trajectory, and zero blood velocity. 38 The dA and A from the equations are measured with the US probe, exerting mechanical pressure on the artery. It is plausible that this influences both arterial diameter and mechanical properties of the artery. Also, dP is measured as the difference between diastolic and systolic pressures as measured at the wrist. This pressure difference is expected to differ from the pressure at the CCA. Also, ρ is based on the literature values only.
Due to the small sample size, it is not possible to draw sound conclusions regarding the reliability of the DIC method. Only a preliminary discussion of the precision, accuracy, and reproducibility of the method will be given. The precision was already discussed earlier in terms of variability. A BA analysis was performed to assess the accuracy. According to such an analysis, DIC estimates PWV on average lower than US, with wide limits of agreement considered from the clinical perspective (see Fig. 10). More research is needed to determine the accuracy of the DIC method in relation to the reference method, given that the cause of the bias is not fully understood. Possible explanations are a flaw in one of the methods, the entirely different natures of the methods, or the fact that the PWV was not Fig. 10 A BA analysis to assess the level of agreement between the PWV acceleration compared to PWV arterial parameters (a) and between PWV dicrotic notch compared to PWV arterial parameters (b). A range of agreement was defined as mean bias (full line) AE2 times SD (dashed line).
simultaneously recorded. Neither a correlation analysis comparing methods nor an intraclass correlation analysis comparing measurement sessions was performed due to the small sample size. Therefore, reproducibility of the method remains to be investigated. Nevertheless, the values obtained from DIC and US are in the same order of magnitude, and they both are in the expected range considering medical background, age, and gender of the test subjects, suggesting that the values are realistic. 7 Additionally, strains were calculated, for illustrative purposes only (see Figs. 8 and 9). Although strains were acquired in all test subjects, data could not be validated with another technique in this experiment. However, as calculation windows were varied, strain patterns remained very similar as assessed with visual inspection. Also, when more heartbeats were recorded in one session, patterns were periodical for separate heartbeats.
The extended spatial and temporal information acquired with DIC opens possibilities for cardiovascular research purposes. Time-resolved, full-field information of the neck region is useful for creating strategies to measure PWV in a noncontact manner. It allows testing for algorithms to detect progress of the wave optically. It allows testing for a robust measurement approach applicable on all types of patients. Moreover, it sheds light on the mechanical behavior of the neck region, and it can be used as the reference for highly realistic FSI models of the CCA. Or, given that reliable FSI models are present, the arterial characteristics can be estimated through reverse engineering, revealing arterial characteristics.
Moreover, DIC can be further developed as a standalone diagnostic technique for PWV detection. One of the benefits of DIC is that the patient can remain fully dressed during measurements. Also, the actual measurement is noncontact; circumventing artifacts due to probe contact in tonometry and US-based applications. The data allow post hoc determination of the arterial trajectory. Additionally, performing a measurement (application of the speckle pattern and recording several heartbeats with DIC) is not likely to take >10 min when executed by a skilled technician.
However, the technique also has its limitations: time-of-flight methods for the PWV detection are typically sensitive to the effect of wave reflections. Also, special attention is required for the preparation of the surface: the contrast between background and speckles has to be as high as possible, and speckles ideally are as small as possible and are evenly distributed in a random fashion, requiring some experience. Moreover, setup of the camera systems, calibration, and especially processing of the data is time consuming making real-time analysis unfeasible at the moment.
Additionally, in order for an object to be studied, visual access from the camera(s) is required over the entire surface, limiting extensive miniaturization of the setup. Also, in order to achieve focus over the entire object, high illumination is required when some reliefs or substantial out-of-plane displacements are at play, especially at the short shutter times. Shadows and blind corners, caused, e.g., by wrinkles, hairs, impurities of the skin, and skin folds will inevitably lead to artifacts. The current measurements were performed in the healthy young males with no facial hair, low BMI, and low PWV. However, it remains to be tested how DIC will perform on skin with skin folds, wrinkles or when there is a prominent presence of facial hair. In order to make the skin surface smoother, a white ground layer capable of filling cracks and pores could be used, and patients could be asked to shave themselves before a consultation. Also, pulsation could be hard to detect when fat content in the neck is high, or when blood pressure is very low. In this case, required sensitivity of the technique could be acquired with a finer speckle pattern. However, the limits of the technique remain to be explored in future work.
Some measures could be made to make the technique more applicable and more user-friendly. The ground layer could be enriched with a fluorescent molecule. By illuminating the patch with a laser and recording with an appropriate filter for the camera, the use of large illumination equipment can be avoided. Or, in order to skip the process of applying two layers of paint, a printed fluorescent pattern could be applied on the skin, e.g., with a rub-on tattoo. The cameras could be triggered by electrocardiogram (ECG) to reduce the amount of images to be analyzed to a minimum. If 3-D information is not required, variants of the technique can be developed with only one camera and even a projected pattern, reducing the cost, reducing the required computational power, and increasing the user friendliness of the setup. In the presence of background movements, principal strains could be used for the estimation of PWV.
It is obvious that the time-resolved, full-field assessment of movement and strains of surfaces creates opportunities for all kinds of biological and nonbiological applications. The technique requires at most two cameras, a bright light source, and the application of some paint, and is physiologically harmless. Moreover, a setup is relatively small and lightweight, and as such very versatile and mobile. There are no interferometric stability requirements, rendering this technique usable for the applications outside the laboratory. Consequently, this technique has the potential to find its way in other biological studies where motion of surfaces is to be assessed as a function of time, for example, study of the movement of eardrums for different frequencies and amplitudes, 39 dynamics of leaf growth, modeling of biting mechanics of stag beetles, and assessing complex biomaterials in general, amongst other possible applications.