Fast segmentation of the femoral arteries from 3 D MR images : A tool for rapid assessment of peripheral arterial disease

Purpose: The peripheral arterial disease is a powerful indicator of coexistent generalized atherosclerosis. As plaques in femoral arteries are diffused and can span a length of 30 cm, a large coverage of the arteries is required to assess the full extent of atherosclerosis. Recent development of 3D black-blood magnetic resonance imaging sequences has allowed fast acquisition of images with an extended longitudinal coverage. Vessel wall volume quantification requires the segmentation of the lumen and outer wall boundaries, and conventional manual planimetry would be too time-consuming to be feasible for analyzing images with such a large coverage. To address this challenge in image analysis, this work introduces an efficient 3D algorithm to segment the lumen and outer wall boundaries for plaque and vessel wall quantification in the femoral artery. Methods: To generate the initial lumen surface, a user identified the location of the lumen centers manually on a set of transverse images with a user-specified interslice distance (ISD). A number of geometric operators were introduced to automatically adjust the initial lumen surface based on pixel intensity and gradient along the boundary and at the center of each transverse slice. The adjusted surface was optimized by a 3D deformable model driven by the local stiffness force and external force based on image gradient. The optimized lumen surface was expanded to obtain the initial outer wall surface, which was subsequently optimized by the 3D deformable model. Results: The algorithm was executed with and without adjustment of the initial lumen surface and for three different selections of ISD: 10, 20, and 30 mm. The segmentation accuracy was improved in a statistically significant way with the introduction of initial lumen surface adjustment, but was insensitive to the ISD setting. When compared with the manual segmentation, the settings with adjustment have, on average, mean absolute differences (MADs) of 0.28 and 0.36 mm, respectively, for lumen and outer wall segmentations, which are significantly lower than those obtained when the adjustment operators were not applied (MAD= 0.43 and 0.59 mm for lumen and outer wall segmentations). The algorithm took about 1% of the time required for manual segmentation to complete segmenting the whole 3D femoral artery. Conclusions: The proposed semiautomatic algorithm generated accurate lumen and outer wall boundaries from 3D black-blood MR images with few user interactions, thereby allowing rapid and streamlined assessment of plaque burden in the femoral arteries. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4916803]


INTRODUCTION
Peripheral arterial disease (PAD), caused by atherosclerosis in the peripheral arterial bed, is a critical worldwide health issue, affecting 27 × 10 6 people in Europe and North America and about 3% in the general Chinese population. 1 The prevalence is significantly elevated for subjects exposed to risk factors such as smoking (15.3%), 2 hypertension (8.7%), 3 and diabetes (16.7%). 4However, PAD is an underdiagnosed and undertreated disease.Late diagnosis of PAD may result in critical limb ischemia, which can lead to nonhealing wounds, gangrene, and eventual limb amputation.More importantly, PAD is an indicator of coexistent generalized atherosclerosis in other vascular beds, which significantly increases the risk of myocardial infraction and stroke for both asymptomatic and symptomatic PAD patients. 5,6Development of imaging tool for detecting PAD would have an enormous impact in reducing the mortality risk of PAD patients, especially symptomatic PAD patients who would likely to benefit from pharmaceutical treatments.
Early investigations of PAD involved imaging the superficial femoral artery (SFA) using x-ray angiography. 7,8How-ever, angiography provides only a two-dimensional (2D) longitudinal view of the lumen.Since no information of the vessel wall can be extracted from x-ray angiograms, angiographic studies are limited to the analysis of stenosis, luminal irregularity, and symmetry. 9Intravascular ultrasound (IVUS) has been shown to be useful in assessing vessel wall thickness, composition, and the degree of luminal narrowing in different segments of the superficial femoral artery.However, both imaging techniques are invasive and not suitable in studies involving asymptomatic subjects and in longitudinal monitoring of PAD patients.
Recent advances in magnetic resonance imaging (MRI) of the vessel wall have allowed direct plaque quantification in the SFA. 10,11The large plaque volume in the femoral artery can be measured reproducibly, leading to decreases in sample sizes and duration required to detect plaque progression/regression with statistical significance.Figure 1 shows an example 3D femoral MR image with manual delineations.As the SFA is long with plaques up to 30 cm in length, 12,13 3D fast black-blood sequences have been developed to attain the longitudinal coverage required to assess SFA atherosclerosis within reasonable scan time. 14The extended coverage afforded by 3D black-blood sequences, however, has posed a challenge in image analysis required for vessel wall and plaque burden quantification.In our experience in generating surrogate ground truth to validate the semiautomated vessel-wall-plus-plaque quantification of the current and a previous study, 15  70-80 min were required to segment arterial lumen and outer wall on 30 axial image slices per artery with 10 mm interslice distance (ISD).The reason of choosing a large interslice distance was to allow the manual segmentation to be completed within a reasonable time frame.However, lumen and outer wall segmentation on axial slices with such a large interslice distance are not enough for reproducible and accurate quantification of vessel wall and plaque burden.Segmenting every axial slices of a 3D MRI acquired in this study, which has a slice thickness of 1 mm, would take at least 700 min.Obviously, a streamlined and validated approach is required for efficient quantification of vessel wall and plaque burden in these 3D MR femoral images with extended coverage.
To the best of our knowledge, only two prior semiautomatic algorithms 15,16 have been described to segment the lumen and outer wall boundaries from 3D black-blood MR images of the femoral artery with extended coverage.Our group has developed a 2D slice-based deformable segmentation approach, 15 in which a user initialized the boundary on the first axial image.The second axial slice was then registered with the first slice using an optical flow registration algorithm. 17he transformation matrix obtained in this registration process was used to transform the boundary on the first image.The transformed contour propagated to the second image and was subsequently deformed and optimized using a B-spline snake model. 18The optimized contour was in turn propagated to the adjacent slice until all axial images were segmented.This method required extensive user interaction to guide the lumen segmentation in regions where blood suppression was poor or where the artery was severely stenotic, thus limiting its efficiency and introducing operator variability.Also, in these problematic regions, longitudinal continuity of the segmented boundaries may not be maintained although registration of adjacent axial slices was performed in the process of contour propagation.In this case, further user interaction on the curved planar reformatted (CPR) view was required to maintain longitudinal continuity.Ukwatta et al. 16 enforced longitudinal continuity of the segmented outer wall and lumen surfaces by maximizing regional overlap between adjacent transverse slices after the stack of contiguous contours was reoriented according to the medial axis of the vessel.The amount of regional overlap alone cannot properly define the relative positions of the two contours if the luminal diameter differs greatly between adjacent slices, such as in stenotic regions depicted in Fig. 2. In this case, the regional overlap is maximized as long as the stenotic contour is within the normal contour.Enforcement of the longitudinal continuity in the original 3D space would be a more robust solution in this situation.However, the spatial information in the original 3D space was removed with the slice-wise reorientation.
In this paper, we describe a segmentation framework that involves much less user interaction and directly optimize the lumen and outer wall surfaces in the original 3D space.This framework requires a user to (a) locate approximate centers of the lumen on a set of transverse images with ISD of 10, 20, or 30 mm and (b) manually segment the lumen and outer wall boundaries on two ends of a 3D image, which are used to estimate how much the optimized lumen surface should be expanded in order to generate the initial outer wall surface.The algorithm does not require any further user interaction after the initial lumen surface has been constructed.Since the 3D deformable model may not be able to drive the initial lumen surface to produce accurate segmentation result if the initial surface is not close enough to the desired boundary, we propose and validate a set of geometric adjustment operators that refine an initial lumen boundary based on the intensity and the image gradient along the boundary and at the lumen center before it is optimized by the 3D deformable model.After the lumen surface is obtained, the initial outer wall surface is constructed by expanding the optimized lumen surface; the outer wall surface is subsequently optimized by the 3D deformable model.The accuracy of the segmented surfaces was evaluated using the same metrics described in Chiu et al. 15 Based on these metrics, we address the two major issues regarding our algorithm: (a) How much is the accuracy of outer wall, lumen, and the vessel-wall-plus-plaque measurements improved after the introduction of the adjustment operators?(b) Are the segmentation results sensitive to the ISD with which lumen centers are manually located?The answer to this question helps us to make a conclusion on the choice of the optimal ISD between manually identified lumen centers that would balance the amount of user interaction and the level of segmentation accuracy.

METHODS
Figure 3 shows the flowchart of the proposed algorithm.The lumen was first segmented in the following three steps: (a) The initial lumen surface was first constructed using B-spline interpolation.The initialization is described in Sec.2.A.(b) A set of adjustment operations were defined to drive the initial surface closer to the desired boundary.These operations are described in Sec.2.B.(c) The adjusted initial surface of lumen was then optimized by a dynamic 3D deformable model, which is described in Sec.2.C.The segmented luminal surface was then expanded and served as the initial surface for the outer wall segmentation, which was subsequently optimized by the dynamic deformable model described in Sec.2.C.

2.A. Initialization of lumen surface
To construct the initial luminal surface, an observer segmented two contours, one on the most proximal and the other on the most distal end of femoral artery within the image volume to be segmented.The user then identified a set of lumen centers on transverse slices of the femoral image.Since the femoral artery is a topological cylinder, the luminal surface can be parameterized by (u,θ), where u is the longitudinal coordinate and θ is the circumferential coordinate representing the angle of surface points with respect to the centerline of the artery.Suppose the arterial surface is represented by p = f (u,θ) and the most proximal and distal contours are equipped with longitudinal coordinates u = 0 and u = L, respectively.Points in the set { f (0,θ), f (L,θ) : θ ∈ [0,2π)} were provided by the observer's initialization of the most proximal and distal contours (see Fig. 4).For the slice with longitudinal coordinate equal to u i that has a lumen center C(u i ) identified by an observer, points on the lumen surface were estimated by the following equation: where r(u i ,θ j ) is the distance between the lumen center on the axial slice with longitudinal coordinate u i [i.e., C(u i )] and the point on the arterial surface with circumferential coordinate θ j .Using Eq. ( 1), eight sample points were obtained at regular intervals of θ for each axial image with manually identified lumen center.B-Spline interpolation was then applied to obtain a dense set of sampled points with resolutions of u and θ equal to 0.4 mm and 0.4 rad, respectively.A surface mesh was then constructed from these sampled points using the Delaunay algorithm implemented in  (The MathWorks, Inc., Natick, MA).

2.B. Lumen contour adjustment
Deformable models require an accurate initial surface, but the initial surface constructed in Sec.2.A may not be accurate enough.In this subsection, we describe a method to improve the accuracy of the initial surface based on the pixel gray level and local geometric structure of the femoral artery on a slice-by-slice basis.We introduce three operators to transform contours on transverse slices: translation, f t , dilation, f d , and contraction, f c , operators, where x is a vertex on the initial contour, c o is the centroid of the initial contour, d is the translation direction vector, τ is a constant controlling the amount of translation, and λ and η are scaling coefficients, with λ > 1 and 0 < η < 1. Figure 5 depicts the effects of these three operators.An automated algorithm was developed to select which operators should be applied to transform initial lumen contour for each transverse slice as described below.Blood inside the lumen appears black in black-blood MRI, and the gray level range of the lumen can be estimated using We defined R v to quantify the percentage of pixels along the initial lumen contour that have gray levels within the range [v lower ,v upper ].The gradient magnitude along the initial lumen contour is expected to be large.We defined R g as the percentage of pixels along the initial lumen contour that have a gradient magnitude higher than a threshold ϵ g .Low R v or R g suggests that a transformation may be required to drive the initial lumen contour closer to the desired boundary.
The gradient magnitude at the centroid, denoted by ∥g o ∥, is expected to be small.Whether or not ∥g o ∥ is smaller than the threshold, ϵ m can be used as a factor in determining the appropriate transformation.In summary, three quantitative metrics are involved in making the decision on which operator should be used in adjusting the position of the initial lumen contour.Two are properties of the boundary, (a) R v and (b) R g , and one is the property of the centroid (c) ∥g o ∥.Let ϵ v , ϵ g , and ϵ m be thresholds for R v , R g , and ∥g o ∥.These thresholds will be optimized by training.We defined the following rules to determine the operator that was applied to the initial contour: (1) No operator was applied if either of the following conditions hold: • R v ≥ ϵ v and R g ≥ ϵ g [Fig.5(a)]; • R g ≥ ϵ g and ∥g o ∥ < ϵ m : In this case, the requirement R v ≥ ϵ v is relaxed, but requirements for the centroid are added.This condition accommodates the situation when the initial contour is poorly aligned with, but close to, the desired boundary.The subsequent 3D deformation is expected to be able to correct this poor alignment [Fig.5(b)].
(2) The translation operator was applied if the following condition holds: • ∥g o ∥ ≥ ϵ m : If the initial contour is shifted with respect to the desire boundary as in the case depicted in Fig. 5(c), the centroid of the initial contour is associated with a high gradient magnitude because it lies near the edge of the desired boundary.In this case, a translation is required.(3) The dilation operator was applied if the following condition holds: • ∥g o ∥ < ϵ m and R v ≥ ϵ v and R g < ϵ g : If the initial contour is inside the desire boundary but too small to capture the whole lumen as depicted in Fig. 5(e), all pixels along the initial boundary fall inside the range [v lower ,v upper ], giving a high R v , but boundary points of the initial contour are associated with low gradient magnitude.The dilation operation is therefore performed to expand the contour.(4) The contraction operator was applied if the following condition holds: • ∥g o ∥ < ϵ m and R v < ϵ v and R g < ϵ g : The initial contour has this property if its centroid is inside the desired boundary, but the initial boundary does not touch the lumen as depicted in Fig. 5(g).The initial contour should be contracted in this case.
Since there are two possibilities associated with each of R v , R g , and ∥g o ∥, there are in total eight possible scenarios.The rules described in the above list are summarized in Table I.
It remained to define the direction vector d of the translation operator.3D initial lumen surface generated by B-spline interpolation tends to be longitudinally smooth and fails to capture the bends of the femoral artery with sufficient accuracy.These locations are where translations are usually required (Fig. 6).We indexed transverse slices in ascending order toward the superior direction.If translation was determined to be needed at Slice i, we calculated the relative position between the contours at Slice i and i − 1 and T I. Choice of adjustment operator based on three parameters extracted from the initial contour: R v , R g , and ∥g o ∥.N, No operation; T, Translation; D, Dilation; C, Contraction.
F. 6.The directional vector, d, associated with the translation operation introduced in Eq. ( 2).d is defined as the sum of p and −g , where p represents the positional difference between the centroids of the boundaries on Slice i, denoted as O i , and Slice i − 1, denoted as O i−1 , and g is the image gradient vector evaluated at O i and projected on Slice i.
denoted it as p, where O i and O i−1 , respectively, represent the centroids of contours at Slices i and i − 1.Since we would like to drive the initial contour toward the lumen, the translation should direct to a region with lower gray level.Thus, we computed the gradient at O i , projected it onto Slice i, and denoted it by g.The translation direction vector, d, was defined as the average direction between p and −g (Fig. 6), The adjustment algorithm was executed iteratively until no operation was required or the number of iterations reached 200.

2.C. Dynamic 3D deformable models
We used 3D deformable models based on the formulation by Terzopoulos and Metaxas. 19,20Suppose the position of a point on the model is x and the displacement is d.The deformable dynamics can be described by the first order Lagrangian mechanics where ḋ = ∂d/∂t, K is a stiffness matrix, and f ext is the external force.K characterizes the elastic properties of the model, which were derived from a local deformation strain energy ε m (Refs.19, 21, and 22) that maintains the continuity of the model surface.ε m can be expressed as where ω 0 weights the displacement magnitude and ω 1 weights the local variations of the displacement.Equivalently, ε m (d) can be expressed as where σ and ε are the stress and strain vectors, respectively.The relationships between ε and d, σ and ε are where in Eq. ( 10), the partial differential operator ∂ was defined as and the symmetric matrix D derived from local deformations is Denoting S as the basis matrix spanning the finite element nodal shape functions, the local displacement d can be approximated in terms of the basis functions as where q d is the vector of local degrees of freedom.Substituting Eq. ( 14) into (10), and then the result into Eq.( 9) yields The stiffness matrix K used in the deformable model is The external forces f ext were defined based on image gradient as follows: where I(x) is the original image, w e is a positive weighting parameter, G σ (x) is a 3D Gaussian function with standard deviation σ, ∇ is the gradient operator, and * is the convolution operator.The model deformation was updated as follows, with the superscript indicating the iteration index, In this paper, deformation terminated when the maximum displacement was less than 0.004 mm or the number of iterations reached 500.

2.D. Laplacian smoothing
We used the Laplacian filter 23 to smooth the optimized surface generated by the 3D deformable model.In the following equation, the position of the vertex p i of the mesh in the nth iteration of the smoothing algorithm was denoted as x n i .The position of p i in the (n + 1)th iteration is given by where p j are vertices with edges connected to p i and λ s is a user-specified weight.According to Eq. ( 22), the new position of p i is its old position translated by the weighted sum of the Euclidean distance between p i and its neighbors.In this study, λ s was set to be 0.01, and the algorithm was executed repeatedly for 100 iterations.

3.A. 3D MERGE acquisition
Six 3D MERGE images acquired for five subjects were used to evaluate the proposed algorithm.Among the subjects, we evaluated the left and right arteries of one subject, while evaluating only one side of the others.
The imaging protocol has been described elsewhere 15,24 and is briefly summarized here.MR images were acquired using a 3 T MRI system (Achieva, Philips Medical System, Best, the Netherlands).The 3D MERGE sequence was implemented using MSDE preparation and spoiled segmented FLASH readout with centric phase encoding.The femoral artery images were acquired using two stations with fieldof-view (FOV) 400 × 40 × 250 mm to cover up to 500 mm longitudinally with isotropic voxel size of 1.0 mm (zerointerpolated to 0.5 mm).Total scan time for one femoral artery was around 7 min.The imaging parameters were TR = 10 ms, TE = 4.8 ms, flip angle = 6 • , turbo factor = 100, and one excitation (NEX).We only analyzed the superficial femoral artery starting from the femoral bifurcation and proceeding to the end of femur.The entire length is approximately 300 mm for each artery.A sample artery with manual segmentation is shown in Fig. 1.

3.B. Evaluation metrics
Distance and area-based metrics were used to evaluate the accuracy of the segmentation.For computing the distancebased metrics, two contours were first matched on a pointby-point basis by a symmetric correspondence technique. 15istance between each pair of points on the two contours was computed and denoted by {d 1 ,d 2 , ...,d N }, where N is the total number of corresponding pairs.Two distance-based metrics T II.Parameters and their optimized values for the adjustment operators.S 1 is the set of optimized parameters obtained in the first training session in which one artery was used for optimization and S 2 is the set of optimized parameters obtained in the second training session in which two arteries were used for optimization.
Area-based metrics were used to compare the area enclosed by the algorithm and the manual segmentations.With the regions enclosed by the algorithm and the manual segmentations denoted by A s and A m , three area-based metrics were defined to summarize the difference between the two contours, which are the percent area overlap (AO), the Dice similarity coefficient (DSC), and the percent area difference (AD), 15,25 where | • | denotes the area of the operand.

3.C. Parameters of the segmentation
The following scheme was used to optimize the parameters for the adjustment operators and the dynamic deformable models.The parameters were initialized empirically, and optimized sequentially by alternating optimization, i.e., changing a single parameter at a time while holding other parameters unchanged.The parameter corresponding to the highest overall AO was chosen as the optimal value.This method is a local optimum search rather than a global one.Thus, the optimized parameter values are not guaranteed to be global optimum.
To evaluate the sensitivity of the parameters to the choice of the arteries used in training, we carried out two training sessions.In the first session, we randomly chose one artery for optimization, obtaining one set of optimized parameters, denoted by S 1 .In the second session, we added another artery randomly to the training set and obtained another set of optimized parameters, denoted by S 2 .Table II shows the parameter values for the adjustment operators.Table III shows the optimized parameter values for the dynamic deformable models.An evaluation on the effect of the use of different parameter sets will be performed in Sec.4.A.

3.D. Validation
The validation consists of two parts: In Sec.3.D.1, the accuracy of the lumen and the outer wall boundaries was evaluated separately based on the distance-and area-based metrics described in Sec.3.B; in Sec.3.D.2, the accuracy of the algorithm in measuring vessel wall and plaque burden was evaluated.Because the vessel wall is in close proximity to neighboring muscle tissues that have similar intensity level to the vessel wall [Figs.1(b) and 1(c)], even an expert observer may have difficulty in segmenting the vessel wall reproducibly.Therefore, the segmentation accuracy of the algorithm should be interpreted with this observer variability taken into account.To quantify this observer variability, an expert observer segmented the lumen and the outer wall on every tenth axial image twice.We randomly chose the boundaries segmented in one session as the gold standard, and T III.Parameters and their optimized values for the 3D deformable model.S 1 and S 2 are two sets of optimized parameters obtained in two training sessions.See Sec.3.C for details.K is the stiffness matrix defined in Eq. ( 16).refer to the segmentation method producing these boundaries as M1 ("manual first session").The boundaries segmented in the second session were denoted by M2 ("manual second session").Metrics described in Sec.3.B were applied to evaluate boundaries generated by M2 and all algorithm segmentation settings against those generated by M1 as described in Sec.3.D.1.By comparing metrics associated with M2 and algorithm segmentation, the accuracy of the algorithm segmentation can be assessed in relation to segmentation variability.

3.D.1. Assessment of lumen and outer wall segmentations
The initial lumen surfaces were constructed as described in Sec.2.A, adjusted as described in Sec.2.B, optimized by the 3D deformable model described in Sec.2.C and smoothened by the Laplacian smoothing filter described in Sec.2.D. Experiments were conducted to investigate the effects of the following three factors on the accuracy of the algorithm segmentation: (a) distance between adjacent manually identified lumen centers, (b) the application of the automatic initial surface adjustment described in Sec.2.B, and (c) the position of the transverse slice relative to the slices where luminal centers were manually identified.To model the effect of these three factors, we defined the following three sets of options that are required to be specified in each algorithm segmentation trial: Distance between adjacent manually identified lumen centers was chosen from 10, 20, or 30 mm.

Option B:
Segmentation was performed with or without the application of the automatic initial luminal surface adjustment."Adj" and "WOAdj" are, respectively, used to denote the options of applying and not applying the adjustment algorithm.

Option C:
Evaluation was either performed on the transverse slices where lumen centers were manually identified [Fig.7(c)] or on slices without manually identified centers [Fig.7(b)].The former option is labeled as "On" and the latter "Btw." The algorithm segmentation was therefore executed in a total of 12 settings (3 possibilities for Option A × 2 possibilities for Option B × 2 possibilities for Option C).A setting was denoted as ⟨Option A⟩-⟨Option B⟩-⟨Option C⟩ where ⟨Option A⟩ = 10, 20, or 30, ⟨Option B⟩ = Adj or WOAdj, and ⟨Option C⟩ = On or Btw.As shown in Fig. 7(a), on each axial slice where manual segmentation was performed, the lumen surface generated in each setting was resliced and compared against the manually segmented contour based on the distanceand area-based metrics described in Sec.3.B.
Lumen contours generated in each of the 12 settings described above were radially expanded to produce the initial contours of the outer wall.The distance that the lumen In (b), evaluation was not performed on slices where lumen centers were manually identified, whereas in (c), evaluation was done on slices with manually identified lumen centers.It is important for a 3D segmentation algorithm to produce accurate segmentation not only on slices with manual center identification but for the whole surface.Comparison between results obtained using the Btw and On settings allows us to assess whether this is the case for our 3D algorithm.contours were expanded was determined by calculating the average wall thickness of the manually segmented outer wall and lumen contours at the two ends of the 3D image.Note that when segmenting the outer wall, the external forces f ext pointed to the brighter region [i.e., the direction of f ext was opposite to that defined in Eq. ( 17)].
Although the adjustment operators were not applied in outer wall segmentation, the initial outer wall boundaries were an expanded version of the corresponding deformed lumen boundaries, and thus adjustment for lumen boundaries had an effect on the outer wall segmentation.The outer wall segmentation generated in each of the 12 settings was also assessed using the distance-and area-based metrics described in Sec.3.B.

3.D.2. Evaluation of vessel wall thickness and area
Rapid and streamlined plaque burden quantification in the femoral artery was the main drive behind developing the proposed segmentation algorithm.Therefore, the accuracy of vessel wall plus plaque burden measurements obtained using different settings should be assessed.Outer wall and lumen surfaces were produced in 12 different settings according to Sec. 3.D.1.Each surface was resliced on axial slice where manual segmentation of outer wall and lumen was performed.Thus, on each axial slice, there are in total of 14 sets of outer wall and lumen boundaries: two sets produced by manual segmentation (M1 and M2) and the remaining 12 sets gener-ated algorithms.Two distance-based metrics and one areabased metric for vessel wall and plaque burden quantification were computed for the 14 settings on a slice-by-slice basis: (a) mean thickness (t mean ), (b) maximum thickness (t max ), and (c) wall area (WA).t mean and t max were defined in the same way as M AD and M AX D except that the correspondence relationship was established between the lumen and outer wall boundaries.WA was the area enclosed by the lumen and outer wall boundaries.Then, the mean thickness difference (∆t mean ), maximum thickness difference (∆t max ), and wall area difference (∆WA) were computed by subtracting the measurement associated with M1 from the corresponding measurement for M2, and for the 12 settings on a slice-by-slice basis.Mean values of these metrics help establish whether there is bias in quantifying vessel wall and plaque burden using the proposed segmentation algorithm.The differences were also assessed by the Bland-Altman plots. 26Another evaluation metric is the root-mean-square error (RMSE), which considered the mean absolute difference instead of the signed difference.

4.A. Effect of segmentation parameters
One artery was chosen and segmented by the 20-Adj-Btw setting of the algorithm twice, once using the S 1 and once using the S 2 parameter set tabulated in Tables II and III.The lumen and outer wall segmentations in a total of 29 transverse slices generated in each of the two segmentation trials were evaluated using the distance-and area-based metrics described in Sec.3.B.The means and standard deviations of the metrics generated by each trial were tabulated in Table IV.T-tests were performed to compare the metrics associated with S 1 and S 2 and no statistically significant difference was found, giving evidence that segmentation performance was not sensitive to whether S 1 or S 2 was used in the algorithm.The experimental T IV.The means and standard deviations of the distance-and area-based metrics associated with parameter settings S 1 and S 2 for lumen and outer wall segmentations.

Lumen
Outer wall results reported in this paper are all generated using the parameter set S 1 obtained using one artery, thereby leaving the remaining five arteries for validation.

4.B. Evaluation of lumen and outer wall segmented boundaries
Tables V and VI show the distance-and area-based metrics associated with M2 and the 12 algorithm segmentations of the lumen and the outer wall boundaries, respectively.Tables VII and VIII show the results of Tukey's multiple comparison tests performed for each of the distance-and area-based metrics for lumen and outer wall segmentations, respectively.In these tables, segmentation methods were ordered in descending order of the respective metric, and parentheses were used to group subsets of methods that are not statistically significant.Figure 8 shows the segmentation results for three sample slices with (a), (b), and (c), respectively, showing examples with good, average, and lower than average segmentation accuracy.The white contours represent the manually segmented lumen and outer wall boundaries and the black contours represent the algorithm segmentation generated using the 10-Adj-Btw setting.Figure 9 shows three examples to demonstrate the T V. Lumen segmentation results: The mean and standard deviation of the distance-and area-based metrics associated with M2 and 12 algorithm segmentations of the lumen boundaries.The gold standard in the evaluation is M1 (i.e., the first manual segmentation).M2 represents a repeated manual segmentation.Each setting of algorithm segmentation is associated with three options, which specifies (a) distances between adjacent manually identified lumen centers used for initial surface construction (10, 20, or 30 mm), (b) whether adjustment operators were applied (Adj or WOAdj), and (c) evaluation method (On or Btw).Please refer to Sec. 3 effects of lumen surface adjustment operators on the accuracy of the deformed contour.The initial contours, the adjusted initial contours, the algorithm segmentations with and without lumen adjustment, and the manually segmented contours were superimposed on the axial images.The initial surface adjustment scheme is an iterative process and the adjustment results were generated by up to a sequence of 200 operations.VIII.Results of multiple comparison tests on outer wall segmentations generated by the 13 nongoldstandard methods.The groups in the second column of the table were ordered in descending order of the respective metrics.The parentheses were used to group subsets of methods with means that are not statistically significant at the 5% significance level.Tables V and VI show that the lumen and the outer wall segmentation accuracy were greatly improved with the introduction of the initial contour adjustment operations for lumen boundaries, and according to the Tukey's multiple comparison tests, the improvements were statistically significant.When adjustment operators were used, (i) the distance between adjacent manually identified lumen centers and (ii) whether the centers were identified on the axial slices where evaluations were performed did not have a statistically significant effect on the accuracy of the lumen and outer wall segmentations as quantified by most of the distance-and area-based metrics.However, there were two exceptions: (i) MAXDs for lumen segmentations generated with the 30-Adj-Btw and 30-Adj-On settings were significantly higher than other settings with lumen boundary adjustment as shown in Table VII and (ii) AD for outer wall segmentation generated with the 30-Adj-Btw setting was significantly higher than other settings with adjustment operators applied.These two exceptions suggest that increased distances between adjacent manually identified lumen centers may have a slight effect on the accuracy on lumen and outer wall segmentation.In contrast, when the adjustment operators on initial lumen surface were not applied, the accuracy of the lumen segmentations was more affected by the distance between manually identified lumen centers.Table VII shows that 30-WOAdj-On and 30-WOAdj-Btw had a significant difference with 10-WOAdj-On and 10-WOAdj-Btw in MAD, MAXD, and AO.F. 9.The effect of the adjustment operators lumen segmentation accuracy.(a), (b), and (c) show cases where the translation, dilation, and contraction operations, respectively, played the most important role among the three operators on improving the initial lumen boundary.The initial contour, the contour after adjustment and the contour after adjustment and deformation are superimposed on each axial image.The deformed contours generated without adjustment operation were also plotted to show the effect of the adjustment operation on the accuracy of the optimized contour.The manually segmented boundaries were also plotted.

4.C. Evaluation of vessel wall thickness and area
The second to fourth columns of Table IX show the difference between the vessel wall and plaque burden measurements associated with 13 nongold standard methods (M2 and 12 settings of the algorithm) and M1, while the fifth to seventh columns of Table IX show the RMSEs of the 13 nongold standard methods compared with M1 as the gold standard.The measurement bias associated with the 12 settings was not different from that of M2.Among algorithm segmentations, RMSEs for all measurements generated with adjustment operators applied were lower than those generated without the adjustment operators.
Figures 12-14 show the Bland-Altman plots of t max , t mean , and WA, comparing 20-Adj-Btw, 20-Adj-On, and M2.It can be observed from Figs. 12 and 13 that the segmented results of 20-Adj-Btw and 20-Adj-On tended to underestimate t max and t mean if the average measurement was larger than a threshold, which was about 2.0 mm for t max , and about 1.6 mm for t mean .But in most cases, average t max and t mean were below the threshold, which explains the positive bias of t max and t mean for these two settings as shown in Table IX.

4.D. Time requirement
All experiments were performed on Visual Studio ++ platform based on VTK (Ref.23) and ITK (Ref.27) libraries, using an Intel(R) Core(TM) i7-3770 CPU@3.40Hz with 8 GB memory.To generate an initial surface for the lumen of one femoral artery with about 300 axial slices took 1-2 min.The adjustment operators took about 30 s to adjust the initial luminal surface.And the dynamic deformable model took about 3-4 min to optimize the lumen and outer wall surfaces.Hence, about 6 min were required in total to segment one artery.In contrast, manual segmentations of the lumen and outer wall boundaries for one femoral artery on every tenth axial slice (i.e., a total of 30 slices/image) took about 60 min, T IX.The mean, standard deviation, and root-mean-square error of the maximum thickness difference (∆t max ), the mean thickness difference (∆t mean ), and wall area difference (∆WA) associated with manual segmentation (M2) and other algorithm segmentation methods (with and without using the adjustment operators) as compared to the gold standard (M1 and it would take a minimum of 600 min to manually segment the whole 3D femoral artery.

DISCUSSION AND CONCLUSION
The development of 3D MERGE technique has allowed fast acquisition of 3D black-blood MR images, which has made the assessment of the full extent of atherosclerosis in the femoral artery feasible.However, quantification of vessel wall and plaque burden in a 3D image of such an extended coverage requires the development of efficient techniques for extracting the outer wall and lumen surfaces.To address this computational challenge, we introduced an efficient segmentation method based on a 3D deformable model in this paper.Deformable models, however, are sensitive to initialization, and our results presented in Secs.4.B and 4.C and the examples displayed in Fig. 10 gave evidence that the 3D deformable model cannot be directly applied to the 3D MERGE MR images due to the difficulties of segmenting the lumen and outer wall boundaries from 3D black-blood MR images.In order for the deformable model to converge to the desired boundaries, we have developed a novel technique involving incremental and iterative translation, dilation, and contraction operations to improve the accuracy of the initial boundary.
We used a set of distance-and area-based metrics to assess the accuracy of lumen and outer wall boundaries in different experimental settings to study the effects of three different factors on the segmentation accuracy.We also assessed how accurate each setting was on quantifying vessel wall and plaque burden by measuring the mean, maximum wall thickness, and area of vessel wall.First, and most importantly, we evaluated whether the introduction of the adjustment operations on the initial lumen boundary has made improvements on segmentation accuracy by executing the algorithm with and without the adjustment operations.In Sec.4.B, Tukey's tests showed that, for all distance-and area-based metrics used, the lumen and outer wall boundaries with adjustment were statistically more accurate than those produced without adjustment of the initial lumen surface.Second, lumen centers were manually identified for a set of axial images with ISD of 10, 20, and 30 mm.The results showed that segmentation accuracy was insensitive to ISD settings when the adjustment operations for initial lumen contour were performed.Only MAXD for lumen segmentation and AD for outer wall segmentation indicated that the ISD setting of 30 mm generated statistically less accurate boundaries compared to the ISD settings of 10 and 20 mm.Thus, to reduce operator interactions in the segmentation pipeline, an ISD of 30 mm can be used with little detrimental effect on the lumen and outer wall segmentation accuracy, although on the conservative side, we recommend an ISD of 20 mm because the additional initialization time required compared to an ISD of 30 mm is less than 1 min/artery.Third, the 3D segmentation method is expected to produce accurate boundaries for the entire image, not just on the axial slices where lumen centers were identified.We evaluated the segmentation accuracy on the axial slice where lumen centers were manually identified, which was compared to that obtained on axial slices where center identification was not performed.Tukey's tests for the lumen and outer wall segmentations for all distance-and area-based metrics (Tables VII and VIII) did not show statistically significant difference between segmentation accuracies on axial slices with and without lumen center identified, giving evidence that accurate segmentation results were produced for the whole artery, instead of only locations where lumen centers were manually identified.
To the best of our knowledge, only two algorithms have been described previously (Chiu et al. 15 and Ukwatta et al. 16 ) to address the computational challenges posed by the requirements of analyzing 3D black-blood MR femoral artery images with extended coverage.Comparing to our previously published 2D-slice propagation method, 15 the proposed 3D method has a number of major advantages.The 2D-slice propagation method required extensive user interactions to steer the algorithm to the correct lumen boundary in regions where the artery was difficult to segment.Further interaction on the CPR view was required to correct for the longitudinal continuity of the vessel.These two steps required a manual observer to monitor the segmentation for a total of up to 7 min for outlining each vessel.In contrast, the proposed 3D algorithm only requires a user to (i) enter lumen centers on a set of axial images with a predefined ISD and (ii) segment the lumen and outer wall boundaries on the ends of the image for vessel wall thickness estimation.The whole initialization process took 1-2 min depending on the ISD used.Hereafter, steps involved including the lumen boundary adjustment, 3D deformation for the lumen surface, generation of the initial outer wall surface, and 3D deformation for the outer wall surface are all automatic.This workflow does not only save time but also reduces operator variability in segmentation, and therefore, vessel wall and plaque quantification.Second, in the previous 2D-slice-based algorithm, although longitudinal continuity was implicitly enforced by optical flow registration of adjacent axial images, segmented contours on adjacent slices could still be very different in severely stenotic regions and regions associated with poor blood suppression.That was the reason why editing was required in the CPR view.In contrast, no user interaction was required in enforcing the longitudinal continuity in the proposed algorithm, as longitudinal continuity in the 3D domain was automatically taken care of by the stiffness matrix in the 3D deformable model and also subsequently by the Laplacian smoothing operation.
Ukwatta et al. 16 introduced indicator functions that label each pixel on each axial image slice as the wall, lumen, and background regions.Then, a coupled continuous min-cut model was applied to minimize a cost function consisting of three terms that (1) match intensity histogram of the wall, lumen, and background regions of each axial image slice with the intensity histogram priors determined by manual segmentation at vessel ends, (2) promote gradients of label indicator functions in regions where the magnitude of the image gradient is high, and (3) promote spatial consistency between slices.However, term (3) was formulated in a way that requires the alignment of the vessel lumen to a cylinder according to a semiautomatically identified medial axis [Fig.2(a)].The spatial constraint was computed not based on the absolute positions of two adjacent slices, but on their relative positions with respect to the medial axis.In regions such as that depicted in Fig. 2 where a stenotic luminal contour is immediately adjacent to a normal lumen contour, the regional overlap quantified by term (3) is maximum when the stenotic contour is inside the normal contour regardless of the relative position of the stenotic contour [Fig.2(b)].Therefore, in this situation, it is not suitable to quantify longitudinal continuity using a regional overlap metric without considering the absolute positions of the two adjacent contours.There are more than ten stenotic regions in our dataset.Figure 10 shows two examples.Figure 11 also shows several stenotic regions along the two example arteries.Thus, there is a need of a more robust metric to quantify the longitudinal continuity in these regions.In comparison, the optimization of point-wise continuity in the original 3D space by the 3D deformable model is a more robust approach and is expected to generate more accurate segmentation results in stenotic regions.Table X compares the segmentation accuracy achieved by the proposed method with those of the previous two algorithms.We used the 20-Adj-Btw setting to represent the proposed algorithm.The proposed algorithm has an accuracy higher than Ukwatta et al. 16 and comparable with Chiu et al. 15 The time required for initializing the proposed algorithm (30-60 s) was less than that required in Ukwatta et al.F. 15.An example demonstrating a limitation of t mean in assessing the accuracy of algorithm segmentations.The setting 10-Adj-Btw was much more accurate in segmenting the lumen and outer wall than the setting 30-Adj-On according to the MAD metric calculated independently for lumen and outer wall and displayed in the bottom of the figure.However, the difference between t mean associated with 30-Adj-On (1.03 mm) and that with the gold standard M1 (1.04 mm) is smaller than the difference between t mean associated with 10-Adj-Btw (1.14 mm) with that of M1.This is due to the fact that 30-Adj-On undersegmented the lumen and outer wall by similar amount pointed to by the white arrows.Comparison among the three labelled t max values, however, showed that contours generated by 10-Adj-Btw were more accurate.
The main drive of the development of a femoral artery segmentation technique was to facilitate a rapid and streamlined process for vessel wall and plaque burden assessment.Therefore, in addition to evaluating the accuracy of the lumen and outer wall boundaries independently, we also assessed the ability of each setting of the algorithm to generate accurate measurements of mean and maximum vessel wall thickness and vessel wall area.However, we observed in our experiments that accurate vessel wall and plaque measurements do not necessarily reflect the accuracy of the segmented lumen and outer wall boundaries.Figure 15 shows a case where we can determine visually that the segmented boundaries obtained using the 10-Adj-Btw setting are much more accurate than those obtained using the 30-Adj-On setting.When the lumen and outer wall boundaries were evaluated independently, we would arrive at the same conclusion.However, the difference between t mean associated with 30-Adj-On (1.03 mm) and that with the gold standard M1 (1.04 mm) is smaller than the difference between t mean associated with 10-Adj-Btw (1.14 mm) with that of M1.The agreement between t mean measured by 30-Adj-On and M1 was due to the fact that 30-Adj-On undersegmented the lumen and outer wall boundaries by a similar amount (white arrows in Fig. 15).This case demonstrates the limitation of the metrics based on vessel wall measurements and suggests that these measurements should be assessed together with the independent evaluation of lumen and outer wall boundaries before arriving at an objective conclusion regarding segmentation accuracy.
The motivation in quantifying vessel wall and plaque volume in the SFA stem from (i) its prognostic implication to future vascular events and deaths 5,6,28 and (ii) the large plaque volume in the femoral artery that can be measured reproducibly, leading to decreases in sample sizes and duration required to detect plaque progression/regression with statistical significance. 10The development of the proposed 3D femoral segmentation algorithm allows for rapid quantification of vessel wall and plaque burden in a 3D black-blood MR image with few user interactions.The time required to segment the SFA with 300 mm coverage was 1% of the time required for manual segmentation.This tool can be used to generate vessel wall and plaque measurements in the SFA efficiently in large cross-sectional and longitudinal studies for assessment of the effectiveness of antiatherosclerotic therapies and evaluation of vascular risk factors.

F. 1 .F. 2 .
An example femoral MR images.(a) shows a 3D femoral MR image with manual segmentation on every tenth transverse slice.ISD of the 3D MR image is 1 mm.Thus, ISD of manual segmentation is 10 mm.(b) and (c) show a sample transverse slice of the 3D femoral MR image with the manual lumen and outer wall segmentations superimposed.Medical Physics, Vol.42, No. 5, May 2015 Schematic of the longitudinal and transverse views of a stenotic region used for illustrating a weakness of the regional overlap metric described in Ukwatta et al. (Ref.16) (a) The longitudinal view of the artery with the medial axis represented by the dotted line.(b)In this transverse view, the distal transverse slice was aligned with the adjacent slice where the stenotic region is located.The regional overlap metric is maximum when the stenotic contour is inside the normal contour regardless of the relative position of the stenotic contour [i.e., the metric does not change even if the stenotic region (black contour) is moved to the left or right (gray dotted-line contours)].

F. 3 .
Block diagram of the proposed segmentation algorithm.Medical Physics, Vol.42, No. 5, May 2015 F. 4. Schematic showing the two required user inputs for initial lumen surface construction: (i) Center of lumen in a set of axial slices with predefined ISD and (ii) manually segmented boundaries on the two ends of the 3D MR image.

F. 5 .
Illustration of the adjustment operators.(a) and (b) show cases for which no operator is required.(c), (e), and (g) show cases for which the translation, dilation, and contraction operators are required, respectively.(d), (f), and (h) show the adjusted contours corresponding to the cases depicted in (c), (e), and (g), respectively.a prior model which can be obtained by training.Most pixels along the initial lumen contour is expected to have gray level within the estimated range, which we denote by [v lower ,v upper ].

F. 7 .
Illustration of Option C of the algorithm segmentation discussed in Sec.3.D.1.Evaluation is always performed on the slice where manually segmented contours are available.(a) shows axial slices with manually segmented contours.(b) and (c), respectively, show the Btw and On settings.

F. 8 .
Segmentation results for three sample axial images.(a), (b), and (c), respectively, show examples with good, average, and lower than average segmentation accuracy.Algorithm and manual segmentations were, respectively, colored in black and white.

F. 10 .
Segmentation results at two example stenotic regions.(a) and (d) show the 3D luminal surfaces segmented by the proposed algorithm using the 20-Adj-Btw setting.(b) and (c) show the segmented lumen and outer wall contours corresponding to the stenotic region shown in (a).(e) and (f) show the segmented lumen and outer wall contours corresponding to the stenotic region shown in (b).The contours segmented without adjustment were generated using the 20-WOAdj-Btw setting of the algorithm.F.11.Manual and segmentations for lumen and outer wall boundaries of two example 3D MR images.(a)-(d) show boundaries segmented for Artery 1, and (e)-(h) show boundaries segmented for Artery 2. (a) and (e): Manually segmented lumen boundaries.(b) and (f): Manually segmented outer wall boundaries.(c) and (g): Luminal surface generated by the proposed algorithm.(d) and (h): Outer wall surface generated by the proposed algorithm.The 20-Adj-Btw setting was used in algorithm segmentation to generate these surfaces.

F. 13 .
Bland-Altman plots of t mean comparing (a) M2, (b) 20-Adj-Btw, and (c) 20-Adj-On with the gold standard M1.Difference values represent t mean of M1 subtracted from that of algorithm segmentation methods.Lines denoting the mean difference and ±1.96 standard deviations are also shown.

F. 14 .
Bland-Altman plots of WA comparing (a) M2, (b) 20-Adj-Btw, and (c) 20-Adj-On with the gold standard M1.Difference values represent WA of M1 subtracted from that of algorithm segmentation methods.Lines denoting the mean difference and ±1.96 standard deviations are also shown.
(98 s), and significantly less than that required in Chiu et al. (up to 7 min).
.D.1 for details.T VI.Outer wall segmentation results: The mean and standard deviation of the distance-and area-based metrics associated with M2 and 12 algorithm segmentations of the outer wall boundaries.The gold standard in the evaluation is M1 (i.e., the first manual segmentation).M2 represents a repeated manual segmentation.Each setting of algorithm segmentation is associated with three options, which specifies (a) distances between adjacent manually identified lumen centers used for initial surface construction (10, 20, or 30 mm), (b) whether adjustment operators were applied (Adj or WOAdj), and (c) evaluation method (On or Btw).Please refer to Sec. 3.D.1 for details.
Medical Physics, Vol.42, No. 5, May 2015 ). t max and t mean are expressed in mm, and WA is expressed in mm 2 .On with the golden standard M1.Difference values represent t max of M1 subtracted from that of algorithm segmentation methods.Lines denoting the mean difference and ±1.96 standard deviations are also shown.
T X.Comparison with two previous algorithms using the same data set.