Semi-automatic method for segmentation of the internal jugular vein in ultrasound movies evaluated at different body postures

Objective: The collapse of the internal jugular vein (IJV) regulates intracranial pressure (ICP) in upright body positions. The cross-section area (CSA) is therefore of interest when studying the effects of postural changes in various neurological diseases. We have developed a semi-automatic segmentation method, which tracks the CSA of the IJV in ultrasound movies, and evaluated its performance in three body positions (supine, 16°, 71°). Approach: The proposed method utilized post-processing image filtering combined with a modified snake active contour algorithm. The ultrasound movies were retrospectively analysed (n = 231, 3s, 28 fps) based on previously collected data from 17 healthy volunteers. The computed CSAs (CA) from the segmentation method were compared to manually segmented CSAs (MA) in two frames per movie. Tracking performance were evaluated by visual inspection. Main results: In the supine position, 100% of the ultrasound movies were tracked successfully, and the mean of CA-MA was −4.4 ± 6.9 mm2 (MA, 88.4 ± 50.5 mm2). The most challenging movies occurred in upright body posture where tracking success rate was 90% and mean of CA-MA was −1.4 ± 2.2 mm2 (MA, 12.0 ± 11.1 mm2). The semi-automatic segmentations took 55 s to perform on average (per movie) compared to manual segmentations which took 50 min. Significance: Segmentations made by the proposed method were comparable to manual segmentations in all tilt-angles, however much faster. Efficient and accurate tracking of the CSA of the IJV, with respect to postural changes, could help furthering our understanding of how IJV-biomechanics relates to regulation of intracranial pressure in different neurological diseases and physiological states.


Introduction
Tracking the cross-sectional area (CSA) of the internal jugular vein (IJV) over time in different body postures is relevant for several neurological diseases as the collapse of these veins has been shown to regulate intracranial pressure (ICP) in upright body positions [1]. ICP varies with body posture [2][3][4][5], yet it is usually studied in the supine position even though humans spend most of their day in an upright position. The IJV collapse could thus be of great significance in the pathophysiology of diseases such as normal pressure glaucoma, idiopathic intracranial hypertension and normal pressure hydrocephalus as little is known about the postural ICP variability in these patients [6,7]. Furthermore, since the IJV is easily accessible with ultrasound, its size and shape have been studied in search for noninvasive assessment methods for different physiological states and indices such as hypovolemia [8] and central venous pressure (CVP) [9,10]. Consequently, being able to track the CSA of the IJV over time, and in different body postures, could provide further insights in how the biomechanics of the IJV relates to various physiological states and neurological diseases.
As the shape and size of the IJV varies related to the cardiac cycle and respiration, measurements of the IJVs should be performed over time. Ultrasound allows for measurements of shape, flow velocity and the CSA of the IJVs. The spatial and temporal resolution of the ultrasound movies makes it possible to study the dynamics of the CSA in real time through the heart cycle. In order to assess all the dynamic effects of the IJV such as variability of pressure and flow resistance, tracking of the CSA over time is crucial. However, performing frame-by-frame measurements manually is highly time consuming, even for small data sets, suggesting the need for automatic segmentation methods. Such an automatic method would have to overcome several obstacles such as: noise interfering with the boundaries of the IJV, large interframe displacements of the vessel walls (due to low frame rate or fast movements of the vessel wall), incomplete boundaries (shadow obscuring the vessel wall), handling very large differences in size and shape (IJVs vary significantly in different body postures and over time) [11][12][13]. To address these issues, previous studies have proposed semi-automatic segmentation algorithms based on different approaches [14][15][16][17][18][19]. In this study we have investigated the feasibility of applying post-processing image filtering techniques combined with an active contour algorithm, based on the snake algorithm proposed by Kass et al [20]. The proposed method was evaluated against manual measurements for 231 ultrasound movies in three different body positions (0°, 16°and 71°).

Materials and method
Subjects and data acquisition We retrospectively analysed ultrasound movies that were acquired in a previous study [11], including 17 healthy volunteers, 10 women and 7 men (age 45± 9 years), all recruited from a newspaper advertisement. The Regional Ethical Review Board in Umeå approved the study (approval no. 2014/223-31), and all participants gave their informed, written consent. Subjects were excluded if they had any neurological, psychiatric or cardiovascular diseases, had blood pressure above 140/90 or if they were using any medication affecting the cardiovascular or central nervous system. The jugular veins were examined using a GE Vivid E9 ultrasound system with a 9 L linear probe (4-8 MHz) (General Electric Healthcare, Chicago, IL, USA) in brightness-mode. Both jugular veins were assessed in 2-3 different neck levels (located 23±2, 27±3 and 29±3 cm from the bottom of the sternum, respectively). Each movie consisted of 89-90 frames, frame rate 28 fps, adding up to approximately three seconds. Tilt-angles 0°, 16°and 71°(i.e. from supine to sitting) were included in this study giving a total of 231 movies.

Proposed method
The segmentation method was written in MATLAB R2017b (The MathWorks, Inc) and run on a HP Intel (R) Core(TM) i7-6700 CPU 3.40 GHz. An open source script (Ritwik Kumar, 2010, https://se. mathworks.com/matlabcentral/fileexchange/28109snakes-active-contourmodels) was used as a foundation that was further developed and modified. Figure 1 illustrates the process in which the method is executed. Figure 1. Once a movie has been uploaded, the user may choose to post-process the images or change any parameter that governs the segmentation algorithm. All images are then filtered with a Gaussian filter. The initial contour vertices are placed manually, using the mouse cursor, on the first frame before initiating the algorithm. The final contour on each frame is used as the initial points for the subsequent frame.

Setup and initiation
The proposed method provided a setup window (figure 2) where the parameters governing the algorithm were set and two optional features selected, an intensity filter (IF) and a balloon force. The IF, when applied, increased the contrast between the edge and the lumen of the vessel by setting all pixel values less than a pre-defined threshold (defined by the user) to zero. The image forces (explained in the next section) that attract the contour were strengthened in this way. The balloon force was defined as a force acting in the normal direction of the contour which prevented the contour from shrinking where the image forces were weak. The setup window also provided a tool for reviewing the movies, frame by frame, to help choosing the appropriate setup and to study the effects of the IF on subsequent frames.
The variety of shape, noise and image quality of the individual movies gave rise to three different setups that were used: (i) active balloon force without IF (default), when large parts of the vessel wall was eliminated using the IF, (ii) active balloon force with IF, when the vessel wall was still continuous even though having applied the IF, (iii) inactive balloon force with IF, when the geometry of the IJV were relatively unchanged even though having eliminated large parts of the vessel wall using the IF (see figure 3). The frequency used for each setup is presented in the results section for each tilt-angle.
A Gaussian filter was applied on the image sequences, regardless of setup, extending the image forces further outside the edge of the vessel wall (further explained in the next section). The user manually placed the initial contour vertices on the first frame close to the true vessel wall (using the mouse cursor) before initiating the segmentation algorithm. Figure 4 shows an example were the IF has been applied prior to initiating the segmentation algorithm.

Segmentation algorithm
From, [20], the snake is given by minimizing the energy, E, of the contour, V, made up by vertices

E int is the internal energy of the contour making it smooth and is defined as
where v s and v ss denotes the first and second order derivative with respect to s, α and β are constants.
Defining the internal energy this way makes the contour continuous and resistant to stretching and bending depending on the values of α and β. Minimizing the first term will pull the vertices closer to one another. Minimizing the second term will prevent the vertices from being arranged as sharp corners making the contour smooth. E ext is the external energy derived from the image properties and is defined as where I is the image intensity, G s is the Gaussian with standard deviation , s w edge and w line are constants greater than zero. When minimizing (1), the first and second terms in (3) provide image forces that make the contour attracted to edges and to low intensity regions respectively. The external forces in equation (3) are local forces that are present only in the immediate neighborhood of the edge. If the edge of the IJV undergoes large inter-frame movements the initial contour vertices on the subsequent frame might be outside the reach of the pulling forces from the edge. To mitigate this problem a Gaussian filter was applied on the images prior to initiating the algorithm. The Gaussian filter smoothens the edge of the IJV and therefore extends its pulling forces further outside the edge making the algorithm more robust.
The Euler-LaGrange equation for minimization of equation (1) A contour made up by n vertices is described in discrete form by the following equation: Approximating with finite differences, the discrete representation of equation ( Equation (6) can be written for all vertices, and for x and y separately, as the following system of equations is the pentadiagonal matrix obtained from the discrete representation of the internal energy in (6). Equation (7) can be solved iteratively in time by replacing the right-hand side with The image forces at time, t, are approximated as equal to the image forces at t 1.
- (6) and adding a weight parameter, , k for the image forces gives

Making the substitution into equation
Moving x t and y t to the left-hand side the following iterative scheme for x t and y t is obtained: The balloon force is defined as an expanding constant force in the normal direction, N v s , ( ( )) at each vertex of the contour. The normal at each vertex, v x y , The implementation of (10) was provided from an open source script online (Dirk-Jan Kroon, 2011, https:// se.mathworks.com/matlabcentral/fileexchange/ 28149-snake-active-contour?s_tid=prof_contriblnk). Conditions for the balloon force were set so that it would only be active when contour points were placed on pixel intensities, I, below a predefined threshold, T. The threshold was set to T 70 = for setup (i) and T 50 = for setups (ii) and (iii). The thresholds were derived from experimenting with different thresholds on a sample of movies.
The final iterative scheme is written separately for each vertex coordinate as where A is the pentadiagonal matrix, x t and y t are the coordinates of each contour point at the iterative step t, f represent the external forces, and , g , k w b are constants. All parameters were set according to what had given the best performance (by visual inspection) for several different types of IJV movies. Parameters were set to: After the final iterative step, the CSA of the contour was computed. The CSA was defined as where P contour is the number of pixels bounded by the contour, P image is the total number of pixels that make up the image, and A image is the area of the image. The contour points were then rearranged equidistantly with the distance, d, along the contour before being used as initial contour vertices for the subsequent frame. Since the internal forces depend on the distance between the vertices as well as the curvature, values of d were conditioned to vary depending on the area of the contour. This made the internal forces less dependent on size and shape. The conditions for d were derived from brief experiments with different conditions on three movies that varied in size and shape. The distance, d, were not allowed to exceed the length of one pixel (0.111 mm) and was also bounded by a lower limit (N lowlimit ) and a higher limit (N highlimit ) with regards to the total number of vertices. Between N lowlimit and N highlimit , d was determined by the size of the CSA. The conditions that determined d are presented in more detail in table 1.

Performance evaluation
The properties of the segmentation method aimed for in this study were (i) tracking the dynamic course of the CSA, confirmed by visual inspection, (ii) estimating the area as close to manual estimates (golden standard) as possible. Satisfying these two requirements would give an adequate representation of the biomechanics of the IJV.

Area estimate evaluation
Reference areas were collected from a previous study [11], where estimates of the maximum and minimum CSA had been made for all movies manually, giving area estimates in two frames per movie. The CSAs computed by the algorithm (CA) were compared to the manual segmented CSAs (MA) of the same frames. CA and MA were both estimated by equation (12).

Tracking evaluation
The segmentation of each movie was inspected visually and categorized as one of two categories, passed or failed. Movies where the contour was significantly misplaced for considerable parts of the movie, where it broke down, or where it was obvious from visual inspection that the contour did not track the edges of the IJV at all, were categorized as failed. However, the computed area value might still be close to the reference area by chance. The area estimates from these movies were excluded from further analysis for this reason. Figure 5 show snapshots from all the twelve failed movies. Ultrasound movies were not excluded with regards to image quality, challenging shapes (such as non-convex shapes) or for any other reasons.

Frame by frame comparison
The proposed method was also evaluated in a frame by frame comparison, on three randomly selected movies (one from each tilt-angle) for which all 90 frames were manually segmented. The DICE-coefficient [21], area variations over time and the total time for the manual segmentations were studied. The DICE-coefficient, S,

Statistics
The CAs and the MAs were presented in a Bland-Altman plot, including maximum and minimum areas together. The correlation between the MAs and the CAs were computed as well. A student's t-test, with the null hypothesis that the mean of the difference, CA-MA, was non-zero, was made for each tilt angel. The significance level was set to p<0.05. Performance with respect to failed and passed movies was presented for each tilt-angle in a separate table.
To study the inter-rater reliability for the manual measurements, a sub-analysis was made on 60 randomly selected movies, one maximum and one minimum in each tilt-angle for 10 subjects. ICC and its 95%-confidence interval were calculated in MATLAB based on a two-way random effect, absolute agreement, single rater/measurement model.

Results
The proposed method was able to successfully track 219 of 231 movies leading to (95%) in the passed-category and 12 of 231 (5%) movies in the failed-category (figure 5) (table 2). The average time consumption for all movies was 55 s (table 2). Figure 6 shows snapshots of movies of different shapes and sizes that were segmented successfully (all three tilt-angles represented).
Tracking of the IJV was increasingly difficult as the tilt-angle increased, which is shown in table 2, where most of the failed movies are seen to occur in tilt-angle 71°, and none in tilt-angle 0°. Also, as table 3 shows, the frequency in which IF was used increased as the tilt-angle increased demonstrating the importance of an appropriate setup when segmenting movies in tiltangles other than supine.
The CAs correlated well with the MAs giving a Pearson coefficient, r=0.99 ( figure 7). The difference, CA-MA, was negative on average for all three tilt-angles, and is seen to decrease as the area of the IJV decreases (see Bland-Altman plot of figure 7 and table 4). The average CA-MA was small compared to MA for all tilt-angles, −4.4±6.9 mm 2 in supine (MA, 88.4±50.5 mm 2 ) and −1.4±2.2 mm 2 in 71°(MA, 12.0±11.1 mm 2 ). Table 4 show the results in each tilt-angle in more detail as well as the results for all tiltangles together.

Frame by frame comparison
The frame by frame manual measurements took approximately 50 min to perform, compared to the computed estimates which took 55 s on average (table 2). Figures 8-10 show the results from tilt-angles 0°, 16°, and 71°respectively. The average DICEcoefficient were 0.95 for the movies in tilt-angle 0°and tilt-angle 16°, and 0.83 in tilt-angle 71°(figures 8-10(b)). Studying the area/frame plot in each of     of each tilt-angle are presented. The average difference between the two raters was −0.8±3.6 mm 2 (MA, 43.7±43.9 mm 2 ) which can be compared to the mean of CA-MA for all tilt-angles together (table 4), which was −3.2±5.3 mm 2 (MA, 50.4±48.7).

Discussion
For performing the challenging task of automatic segmentation of IJV in ultrasound movies, acquired both from supine and sitting posture, a combination of a high-pass filtering and the snake algorithm proved to be a successful approach. The intensity filter, IF, eliminated noise and prevented the contour from converging to incorrect boundaries created by noise inside or outside the true boundaries of the IJV. The IF extracted distinct boundaries of the vein which generated strong image forces that enabled effective tracking of a great variety of different shapes and sizes of the CSA in all tilt-angles. Robustness was supported by a successful tracking in 95 % of the ultrasound movies. The proposed method was much faster than manual segmentations and the CAs were comparable to MAs in all tilt-angles (table 4). However, the relative error of the area estimates is greater for small CSAs than larger ones. The small but significant underestimation of the MAs (table 4) can be partly explained by the algorithm in which the balance of external and internal forces occurs in the immediate neighbourhood of the inner side of the vessel wall. The underestimation is reinforced further by the intensity filter which creates a sharp edge adjacent to the true edge on the inside of the vessel wall. Segmentations made in the supine position, where the IJV is conventionally assessed in both clinical practice and clinical studies, produced good representations of the IJV biomechanics using the proposed method. The right IJV is closely connected to the right atrium of the heart making it an interesting study object for blood volume and CVP assessments. A characteristic shape of the jugular vein has been associated with hypovolemia [8], and the end expiratory diameter of the IJV has been found to be correlated to CVP [9]. In Nakamura et al [22] they used a semiautomatic segmentation algorithm, which tracked the CSA of the IJV in supine posture and found different IJV indices that correlated with stroke volume and estimated CVP (compression sonography of the cephalic vein [23]). The results of this study support the applicability of the proposed method in these fields of research as well. Studying the shape, size and dynamics of the CSA and their association to different physiological states would require movies recorded over several seconds to be reliable (local extreme values vary considerably in area/frame-plots of figures 8-10(c)). Having access to an automatic tracking tool is therefore essential in the data analysis.
In tilted body postures (16°and 71°), the CSAs were small which made it more challenging to locate the boundaries of the IJV. This was especially true when noise was present or when ultrasound movies had poor image quality. However, choosing the appropriate setup, where IF proved to be crucial  (table 3), it was possible to achieve efficient tracking and CSA estimates close to manual segmentations in most of the ultrasound movies in the tilted body postures. When going from the supine to an upright body posture, the IJV collapses as the hydrostatic pressure drops below the external pressure [12]. This collapse have been shown to take part in the postural regulation of ICP [1], suggesting that further investigation of IJV-biomechanics in upright might increase our understanding of ICP variability in this position. Some caution should be taken when calculating flow resistance, which are inversely proportional to area squared, for small CSAs where minor errors can give significant differences in resistance. But in general, the proposed method could be utilized to further investigate whether relative differences exists regarding, CSA, circumference, tilt-angle of collapse, and pulsatile movements between healthy and disease affected subjects. An experienced user might get better results than a novice user in higher tilt-angles, since it requires more fine-tuned setup and careful review of the movies to achieve acceptable results. Total time consumption per movie might therefore increase as well since it is more frequently necessary to terminate the segmentation prematurely, refine the setup, and redo the segmentation for difficult movies. Common among the difficult, and failed movies, were combinations of poor image quality, obscured vessel walls, rapidly   moving vessel walls and CSAs of very small size, all common features in higher tilt-angles. Table 6 shows previous studies that have presented algorithms for IJV segmentation in ultrasound movies. They primarily present performance with respect to the supine position [14][15][16]19]. While Karami et al [17,18] did measure on several tilt-angles they presented a joint performance and for categories such as shape and image quality. None of the referenced studies present performance with respect to different tilt-angles or the mean differences, which has been the main focus of the current study.
The algorithms in table 6 produced the best results in the supine posture where the IJV is distended and usually has an oval shape (in Karami et al [17,18] they did not show performance in the supine posture exclusively however they presented performance with respect to oval shape [18]). This is true for the proposed segmentation method as well (table 4, figure 8(b)). Even though the studies in table 6 uses different indices for evaluation such as DICE [14,15,17,18], Jaccard index [19], accuracy index [16], the strong results in table 4 support that the proposed segmentation method is comparable to these algorithms in the supine posture. In tilted body postures a comparison is difficult to make since our study is the first to present performance for different tiltangles separately.
The positive results of the segmentation method proposed in this study encourage further development. This should include further optimization of the parameters that govern the active contour. Experimenting with considerably lower values of α and β (that govern the internal energy) may enable having less vertices that build up the contour, thus speeding up the segmentation further by having to solve systems of equations of smaller size. Also, adding a feature where the user may interact with the contour manually with the mouse cursor or with the keyboard, aiding in difficult passages of challenging movies, could reduce the need of terminating segmentations prematurely and increase the accuracy. The MATLAB code from the current study is available from the corresponding author.

Conclusion
In this study we have developed a segmentation method that tracks the vessel wall movements of the IJV, and estimates its CSA from ultrasound movies. The proposed segmentation method was evaluated on ultrasound movies recorded in three different body postures, from supine to upright. By combining postprocessing image filtering with a modified snake active contour algorithm the proposed method successfully segmented a great variety of sizes and shapes of the IJV in all three body postures. It may be useful in various fields of research, such as investigation of postural effects in neurological diseases, or how IJV biomechanics are associated with different physiological states.  [16] ImageJ+cross-correlation 3 3 Supine Ikhsan et al [19] Star-based 16 58 Not specified Karami et al [15] Region growth+active contour 2 2 Not specified Karami et al [17] Lucas-Kanade 14 70 0°-90°K arami et al [18] Adaptive polar active contour 13 65 0°-90°C urrent study Intensity filter+Snake 17 231 0°-71°1