Vector Field Convolution-Based B-Spline Deformation Model for 3 D Segmentation of Cartilage in MRI

In this paper, a novel 3D vector field convolution (VFC)-based B-spline deformation model is proposed for accurate and robust cartilage segmentation. Firstly, the anisotropic diffusion method is utilized for noise reduction, and the Sinc interpolation method is employed for resampling. Then, to extract the rough cartilage, features derived from the real symmetric Hessian matrix are chosen to enhance the cartilage, followed by binarizing the images via an optimal thresholding method. Finally, the proposed VFC-based B-spline deformation model is used to refine the rough segmentation. In the experiments, the proposed method was evaluated and demonstrated on 46 magnetic resonance images (MRI) (including 20 hip joints and 26 knee joints), and the results were compared with three state-of-the-art cartilage segmentation methods. Both qualitative and quantitative segmentation results indicate that the proposed method can be deployed for accurate and robust cartilage segmentation. Furthermore, from the segmentation results, patient-specific 3D models of the patient’s anatomy can be derived, which then can be utilized in a wide range of clinical applications, such as 3D visualization for surgical planning and guidance.


Introduction
Segmentation of the articular cartilage is regarded as an essential and primary step in osteoarthritis (OA) quantitative analysis, and can also help clinicians assess the progress of OA treatment [1].However, articular cartilage segmentation from magnetic resonance images (MRI) is still a challenging task due to the thinness of the cartilage, low contrast to surrounding tissues, and the tiny gaps between two cartilage elements of the same articular structure [2].Consequently, conventional segmentation techniques are often insufficient to achieve accurate segmentation of articular cartilage from MRI datasets.
Thresholding and region growing approaches were first proposed for cartilage segmentation because of their fast processing speed.Hardy et al. [4] proposed a combination of thresholding and manual editing for calculating and displaying the distribution of articular cartilage thickness from MR images.Pakin et al. [5] utilized a mixture of region-growing, labeling, and deformable model-type algorithms for cartilage extraction.
Edge detection methods are also efficient approaches for dealing with cartilage segmentation problems.Mahmood et al. [6] applied various edge detection methods for knee joint articular cartilage segmentation and visualization.Jaremko et al. [7] combined manual segmentation with Canny edge detection to separate the tibial cartilage images of different periods and calculate their cartilage volumes, and then registered the tibial cartilage of different periods to quantify the volumetric changes of knee joint cartilage.Carballido-Gamio et al. [8] used edge detection to segment cartilage, meanwhile, an affine registration technique was proposed to compare cartilage thickness.Ghosh et al. [9] and Grau et al. [10] both used the watershed algorithm for the segmentation of the knee cartilage.Folkesson et al. [11] achieved the automatic segmentation of cartilage and non-cartilage regions using a two-step k-nearest neighbor (k-NN) voxel classifier.Although having the advantage of fast processing speed, most approaches in this category have limitations in that they produce noisy images, especially when the noise is present in cartilage layers of highly congruent joints.
The second category is model-based approach.This type of methods make good use of prior knowledge about the target to overcome the limitations of the first category.The methods in this category can be further divided into two types: statistical-based models [12,13] and deformation models [14].
With relative higher accuracy, statistical-based models have been paid much attention by many researchers in recent years [12].However, these approaches tend to be computation intensive.Lee et al. [15] developed a fully automated method for cartilage segmentation using a hybrid multi-atlas segmentation (MAS) and graph-cut-based region adjustment scheme.Glocker et al. [16] performed the patellar cartilage segmentation by fitting a statistical atlas with MRI data in a non-rigid registration protocol.Solloway et al. [17] utilized an active shape model (ASM) [18] to perform the femoral cartilage segmentation task.Ye et al. [19] presented an improved 3D shape context combined with ASM for accurate extraction of bones and cartilage in MR images.Yin et al. [20] utilized the Layered Optimal Graph Image Segmentation of Multiple Objects and Surfaces (LOGISMOS) framework to simultaneously achieve segmentation of all bones and cartilage surfaces in MR images.Lee et al. [21] proposed a hierarchical scheme for the automatic segmentation of knee joints in MR images, which includes four parts: modified branch-and-mincut algorithm-based bone segmentation, bone-cartilage-interface (BCI) detection, optimization of local shape and appearance probabilities-based cartilage segmentation.
Deformation models, especially for active contour models, have received much attention owing to their high processing accuracy [14].These models and their improved methods have also achieved good performance on cartilage segmentation in 2D medical images.However, the main disadvantage of these methods is that 3D contextual information is not fully utilized because of the absence of 3D features in a slice-by-slice fashion [22].Consequently, the accuracy of the segmentation results may not be satisfactory.Stammberger et al. [23] proposed a B-spline snake model for cartilage thickness measurements and interobserver reproducibility assessment.Lynch et al. [24] combined expert knowledge with the active contour model to segment knee cartilage.Kauffmann et al. [25] segmented cartilage using a semi-automatic active contour model while utilizing a point-based rigid-body registration technique to track changes in cartilage thickness and volume.Fripp et al. [26] presented a hierarchical segmentation scheme that includes automatic bone segmentation using 3D ASM, BCI extraction, and cartilage segmentation from BCI via a hybrid deformable model approach.
Among various active contour model-based methods, vector field convolution (VFC) [27] based approaches have shown good prospects for reliable and robust cartilage segmentation.However, currently, these methods still cannot accurately segment the cartilage from MR images for pathological joints.In addition, other challenging problems of cartilage segmentation have further exacerbated this difficulty, including the thinness of the cartilage, low contrast in comparison to surrounding tissue, and the difficulties to accurately locate cartilage interfaces.
In view of the above problems, a 3D VFC-based B-spline surface deformation model is proposed for accurate and robust cartilage segmentation from MR images in this paper.The main contributions of this study are summarized as follows: (1) a Hessian matrix-based cartilage enhancement method is introduced to improve model initialization; (2) The VFC external force field is extended to 3D for the B-spline surface deformation model to increase segmentation robustness; (3) the B-spline surface deformation model without explicit internal force is presented to accelerate the deformation process and reduce the computational cost; (4) The proposed method is successfully applied to extract cartilage surfaces from 46 MR images.
In the first section, a brief introduction of articular cartilage segmentation is described.In the second section, we depict the workflow of the proposed segmentation framework one by one.In the final section, experimental results are demonstrated and discussed.

Methods
In this section, the details of the proposed method will be depicted step by step until the final results are achieved.The flowchart of the segmentation scheme is shown in Figure 1.

Image Preprocessing
In order to reduce the influence caused by image noise and anisotropy, image filtering and resampling operations were used for this purpose.Although a Gaussian filter is widely used for noise reducing, it may blur the edge of the image and to some extent lose the key edge information.Therefore, the anisotropic diffusion method [28], which is a non-linear filter, was employed and implemented by solving the following nonlinear partial differential equations: where div and ∇ denote the divergence and gradient operator, respectively, and g() is the diffusion coefficient.
The typical raw MR images obtained from hospital are usually anisotropic, meaning that the pixel pitch along the z-axis (∆z) is much larger than the pixel pitch (∆xy) in the projective plane.Thus, the images need to be resampled to make the two values (∆xy and ∆z) as close as possible.Motivated by the fact that Sinc interpolation tends to result in minimum aliasing artifacts for uniform sampling, the Sinc interpolation algorithm [29] is specifically adopted for image resampling in this paper.

Cartilage Enhancement
Before cartilage segmentation, we performed cartilage enhancement to improve its delineation accuracy.Inspired by the classic Frangi vessel enhancement filter [30], we adapted it to form a cartilage enhancement filter via the eigenvalue analysis of the Hessian matrix at multiple Gaussian scales.For a 3D MR image I 0 to be processed, the Hessian matrix at voxel x 0 in scale s is a square matrix made up of all the second partial derivatives of the Gaussian-filtered image I: For the Hessian matrix H 0,s , the absolute values of its three eigenvalues can be sorted in an ascending order, i.e., λ 1 ≤ λ 2 ≤ λ 3 , and the corresponding eigenvectors are V 1 , V 2 , and V 3 , respectively.
It has been shown that the three eigenvalues and eigenvectors of the Hessian matrix H 0,s exactly characterize the key features of the 3D surface [30]: 1.The eigenvalue λ 3 with the largest absolute value and the corresponding eigenvector V 3 represent the intensity and direction of the largest curvature at x 0 , respectively.2. The eigenvalue λ 1 with the smallest absolute value and the corresponding eigenvector V 1 represent the intensity and direction of the smallest curvature at x 0 , respectively.
Therefore, because a voxel belongs to a certain morphological structure, some conditions on eigenvalues of the Hessian matrix must be met.Table 1 lists the conditions for all typical morphological structures in 3D images.Note that the sign of the eigenvalues is closely related to the contrast between the morphological structure and image background.In our study, bright structures on a darker background are assumed, resulting in a negative semi-definite Hessian matrix.Therefore, all the eigenvalues of high absolute values are negative.As shown in Table 1, if a voxel belongs to a plate structure, e.g., cartilage, its eigenvalues should satisfy the following conditions: As in the classic Frangi vessel enhancement filter [30], the eigenvalues of the Hessian matrix are combined to form a plateness measure for cartilage enhancement.The mathematical description is as follows: where are two specially designed geometric ratios as in Reference [30], 3 is a measure of second order structureness [30], α, β, c are parameters for the three ratios, and V 0 (s) is the plateness value in scale s.
To account for the varying thickness of the cartilage layers, a multi-scale enhancement approach was adopted.Specifically, the final plateness value of voxel x 0 was set as the maximum response among all the scales ranging from s min to s max :

Image Binarization
From the previous step, the enhancement of the cartilage was accomplished.In this step, a faster version of the Otsu's optimal thresholding method [31] was employed for image binarization.In this method, a modified between-class variance is proposed, which greatly improves the computational efficiency for searching for the optimal thresholds.In order to separate different types of cartilage, the optimal thresholding was followed by the application of a binary median filter, connected-component labeling with seeds, and hole filling using a closing operator.This image binarization step results in the initial locations of cartilage, from which the VFC-based B-Spline deformation model can be automatically initialized.
Figure 2 shows the rough segmentation results of hip cartilage in an MR image.After the preprocessing and cartilage enhancement step, the voxels belonging to the hip cartilage are identified (Figure 2b).Then, after the image binarization step, the femoral head cartilage (Figure 2c) and the acetabular cartilage (Figure 2d) are separately derived.
The resulting cartilage edges are then assigned as the initial contours for the B-spline surface deformation model in the following steps.

Cartilage Refinement Based on B-Spline Deformation Model
The construction of the 3D B-spline surface deformation model is a critical step for the refinement of rough cartilage segmentation, the whole procedure will be described in detail in the following sections.

External Force Field of B-Spline Deformation Model
After the 3D edge image of the cartilage was roughly extracted, it was used to calculate the external force field, to which the deformation surface model was attracted.In this study, the VFC field proposed by Li [27] was employed as the external force field.VFC greatly expands the influence scope of the gradient vector field, which yields a good attraction to the deformation model for space points farther away, and can also let the model converge on the concave boundaries.In addition, compared with the classical gradient vector flow (GVF) external force [32], VFC fields are more robust to noise and initialization while requiring much less computational cost [27].In this paper, unlike the traditional use of the VFC field of the image for 2D contour extraction, we generalize it for the 3D B-spline surface deformation model as the external force field.The VFC field of a 3D input image I(x, y, z) ∈ R 3 is defined as a vector field v(x, y, z) = [u(x, y, z), v(x, y, z), w(x, y, z)], which is computed as the convolution of the edge map with a vector field kernel k(x, y, z): where ⊗ denotes linear convolution, f (x, y, z) is the 3D edge map of the input image I(x, y, z).
The vector field kernel k(x, y, z) is defined as follows: where m(x, y, z) is the magnitude of the vector field kernel, and n(x, y, z) is the unit vector from (x, y, z) to the kernel origin calculated as: In addition, m(x, y, z) is defined as a decreasing function of the distance from the kernel origin [27]: where γ > 0 controls the decrease, and ε > 0 is a small constant to prevent division by zero.
The 3D VFC field v(x, y, z) can be obtained by convolving the edge map f (x, y, z) with each component of the vector field kernel k(x, y, z) as in Reference [27].For 3D edge images, while initializing the edge model (surface), its VFC field v tends to right towards the edge, resulting in larger values in the immediate vicinity of the edges.

The Definition of B-Spline Surface Representation
In this study, the bi-cubic B-spline surface was adopted for geometric representation of the deformation surface model, as it offers a good tradeoff between shape complexity and computational cost.
Bi-cubic B-spline surfaces can be defined by using B-spline basis functions as in Reference [33].Firstly, one point on the (i, j)-th B-spline surface patch can be defined as follows: where u = [u 3 u 2 u 1] and w = [w 3 w 2 w 1] are parameter vectors, B R ij is the 4 × 4 control point matrix with 16 control points ranging from When parameters u and w vary from 0 to 1, all points on the B-spline surface can be calculated from the above equation.
From Equation ( 9), the (i, j)-th B-spline surface patch can be calculated as follows: where , n denotes the number of surface points in both directions while u and w vary from 0 to 1, and P i,j is the n × n surface patch matrix with surface points ranging from P ij (u 1 , w 1 ) to P ij (u n , w n ).

Deformation Process of B-Spline Deformation Model
In the B-spline surface deformation model, the shape of each spline patch is determined by the positions of its 16 control points.At each step of the deformation process, the control point matrix B R ij is updated, and the B-spline surface then continually converges toward the boundary of the target object.However, the external force field does not function directly on the control points but surface points.Therefore, the positions of the control points cannot be updated directly by the influence of the external force field.Nevertheless, by exploiting the relationship between the control points and surface points, the effects of the external force field on the control points can be indirectly derived from the external force acting on the surface points (as shown in Figure 3). Fi,j+1 Fi-1,j+2 Fi,j+2 Fi+1,j+1 Fi+1,j+2 Fi+2,j+1 Fi+2,j+2 Qi-1,j-1 Qi-1,j Qi-1,j+1 Qi-1,j+2 Qi,j-1 Qi,j Qi,j+1 Qi,j+2 Qi+1,j-1 Qi+1,j Qi+1,j+1 Qi+1,j+2 Qi+2,j-1 Qi+2,j Qi+2,j+1 Qi+2,j+2 As can be derived from Equation ( 9), a B-spline surface patch is determined by its 16 control points.If the points on the surface patch are known, then the 16 control points can be obtained by inverting Equation (10).Specifically, from n × n surface points, the corresponding control points B R ij can be calculated as follows: Therefore, the relationship between the control points and surface points in Equation ( 11) can be utilized to obtain the effect of the 3D external force on the control points.
Assuming that the external force acting on the point P ij (u, w) of the surface is FP ij (u, w), then the external force that exerts on the (i, j)-th surface patch can be represented by an n × n matrix FP i,j , with elements ranging from FP ij (u 1 , w 1 ) to FP ij (u n , w n ).
Assuming that ∆B R ij denotes the changes of control point matrix in each deformation step, then from Equation (11), ∆B R ij can be obtained as follows: where η denotes the step size.Then, during each step of the B-spline surface deformation, the control matrix B R i is iteratively updated as follows: Considering that the overall shape of the cartilage surface is simple and smooth, we employ a simple deformation process for the B-spline surface model as suggested in Reference [33].Specifically, multiple surface patches are first independently initialized, and overlapping portions may exist between different surface patches.Then, under the influence of the VFC external forces, each surface patch is independently deformed based on its control points via Equations ( 10)-( 13) in an iterative fashion.After the deformation process of the B-spline surface model terminates, the final cartilage surface can be obtained from resulting surface points.

Image Data and Parameters
A set of 46 MR images were used in this study, including 20 hip and 26 knee images.The experimental data was obtained from a 3T MR machine (GE Medical Systems, Milwaukee, WI, USA), with the slice plane resolution of 256 × 256 and the in-plane pixel pitch of 0.625 mm, the distance between slices was 1.5 mm, and the number of slices ranges from 106 to 217.The experimental data includes four typical types: normal bone and joint, suspected OA, mild OA, and moderate OA.Eleven of the 26 knee joints were normal, and the knee MR images consisted of two bones (femur and tibia) and two corresponding cartilage elements (femoral cartilage and tibial cartilage), as shown in Figure 4. To obtain the golden standard, two radiological experts manually segmented each joint of all the MR images in a slice-by-slice fashion via the ITK-SNAP toolbox [34], resulting in binary mask images of the target cartilage, which are considered as the golden standard in this study.The experimental program runs in MATLAB, with the computer configuration as follows: Intel (R) Core (TM) i7-4770 3.4 GHz CPU, and 16 GB RAM.

Femur Tibia
The parameters used in our method are set as follows: α, β, c in Equation ( 4) are set to 0.5, 0.5, and half the maximum value of S, respectively, and the scales range from s min = 0.5 to s max = 4.0; γ = 1.7 and ε = 10 −8 for the VFC external force field; η in Equation ( 12) is set to 0.5, and the number of B-spline surface patches used is 120.

Comparison of Different Surface Models
Regarding the extraction of the target surface of 3D images, the traditional method is the deformed surface model algorithm [22].During the iterative evolution, the model needs to calculate the internal force for the deformed surface, so as to ensure that the deformed surface is smooth and would not break or crack.While for the B-spline surface deformation model, the strong implicit constraint on the internal force of the surface makes it unnecessary to separately calculate the internal force for the B-spline surface in the whole evolution procedure, which saves much of the calculating cost.
On the basis of the knee images used in this study, the knee joint cartilage was extracted by both the traditional surface model algorithm and our B-spline surface deformation model, and the accuracy and run time were compared.In order to make a fair comparison, the traditional surface model was also implemented in MATLAB and run in the same computer system.The metric of Dice similarity coefficient (DSC) [35] was used for quantitative comparison in this study, which is defined as the volumetric overlap between segmentation results A and golden standard B: DSC = 2A∩B A+B .Table 2 gives the accuracy comparison between our model and the traditional model.Quantitatively, the proposed model outperforms the traditional one when segmenting tibial cartilage and femoral cartilage.In addition, Table 3 shows comparisons of both the number of iterations and running time when segmenting test images with different number of slices.Compared with the traditional model, our B-spline surface deformation model requires fewer iterations to converge, thus resulting in a much lower running time.In addition, one segmentation example of a hip joint with mild OA is shown in Figure 6. Figure 6b shows the 2D cartilage segmentation results.In addition, another segmentation result for a hip joint with moderate OA is shown in Figure 7. Similarly, Figure 7b shows the 2D cartilage segmentation results.Table 4 shows a comparison of the proposed method with other three state-of-the-art methods for automatic segmentation of knee cartilage.It can be seen from the table that our proposed method outperforms the other three methods on segmenting the femoral cartilage, and achieves a slightly better accuracy with regards to the tibial segmentation.The main reason for this is that it is difficult for our method to accurately segment the regions of high curvature without under-segmentation.Therefore, other methods achieves higher specificity through slight over-segmentation, while our method tends to be under-segmentation.
Both the qualitative and quantitative segmentation results presented above indicate that the proposed method is able to achieve high accuracy in cartilage segmentation.Therefore, with high accuracy and efficiency, the proposed method can be deployed for accurate and robust cartilage segmentation.

Discussion and Conclusions
Cartilage segmentation, that is, to determine the boundary of two adjacent thin structures around the joint, is often complicated by the limited resolution of the image data and pathological cases.In this paper, we have proposed a new VFC-based B-spline surface deformation model for accurate and robust 3D cartilage segmentation from MR images to overcome the deficiencies.Our method consists of three main parts: image preprocessing, rough segmentation of cartilage, and refinement of cartilage segmentation.To verify the effectiveness and clinical applicability of the proposed method, extensive evaluations on 46 MR images were conducted in the experiment.The cartilage obtained by our method covers all the cartilage in the MR images, and our method can segment the cartilage very accurately.
The main advantage of the proposed B-spline surface deformation model is its simplicity and efficiency.In addition, the B-spline surfaces facilitate the calculation of different differential attributes, e.g., surface normal and curvature.Generally, in the traditional deformation models, the internal force is introduced mainly to keep the deformed surface continuous and smooth, especially in regions of high curvature.While in the B-spline surface deformation model, the internal force is no longer necessary for the deformation process, because the continuity and smoothness of the B-spline surface are implicitly constrained by the B-spline formulation.Consequently, our proposed B-spline deformation model without explicit internal force is much simpler than the traditional deformation models, and requires significantly less calculation time.In addition, the Hessian matrix-based cartilage enhancement method greatly improves the accuracy of the following model initialization.Furthermore, the VFC external force field for the B-spline deformation model significantly increases segmentation robustness.
After comparison with the three state-of-the-art methods, our proposed method outperformed the other methods in segmenting the femoral cartilage; meanwhile, a slightly better accuracy was achieved on the tibial segmentation.It should be noted that the compared methods are all based on some prior knowledge, the attainable segmentation accuracy thus largely depends on the training set as well as the appropriateness of feature selection.In contrast, our method is not constrained by the training set when performing cartilage segmentation, and thus can be much more effective when only limited data are available.Li, et al. [36] proposed a region-based level set method for MRI segmentation in the presence of intensity inhomogeneities, which is also designed to simultaneously estimate the bias field for bias correction.However, only the multiplicative component was taken into account in their method, which tends to misclassify adjacent thin structures.Furthermore, it is well known that solving level-set PDEs imposes a significant computational burden.Therefore, when applied to cartilage segmentation, the method may lead to a less satisfactory performance.
Despite having achieved promising results, our method still has certain limitations and needs further evaluation on larger image datasets.When performing cartilage segmentation on joints with severe OA, it is difficult for our method to obtain satisfactory results.Therefore, our follow-up study will focus on accurate and robust segmentation and measurement of cartilage with severe OA.

Figure 1 .
Figure 1.Flowchart of the proposed scheme on cartilage segmentation.

Figure 2 .
Figure 2. Initial segmentation of femoral cartilage and acetabular cartilage from an MR image of hip joint.(a) MRI of a hip joint.(b) Articular cartilage after enhancement.(c) Femoral cartilage after image binarization.(d) Acetabular cartilage after image binarization.

Figure 3 .
Figure 3. External forces acting upon B-spline surface points and control points.

Figure 4 .
Figure 4.An illustration of the knee MR imaging.One 2D slice from sagittal MR sequence with anatomical structures annotated.

Figures 5 -
Figures5-7show three typical cartilage segmentation results for the proposed method.Figure5shows the cartilage segmentation results for a knee joint with suspected OA.The original MRI of the knee joint and the 3D visualization of the femoral cartilage obtained by the proposed method are shown in Figure5a,b, respectively.In addition, one segmentation example of a hip joint with mild OA is shown in Figure6.Figure6bshows the 2D cartilage segmentation results.In addition, another segmentation result for a hip joint with moderate OA is shown in Figure7.Similarly, Figure7bshows the 2D cartilage segmentation results.
Figures5-7show three typical cartilage segmentation results for the proposed method.Figure5shows the cartilage segmentation results for a knee joint with suspected OA.The original MRI of the knee joint and the 3D visualization of the femoral cartilage obtained by the proposed method are shown in Figure5a,b, respectively.In addition, one segmentation example of a hip joint with mild OA is shown in Figure6.Figure6bshows the 2D cartilage segmentation results.In addition, another segmentation result for a hip joint with moderate OA is shown in Figure7.Similarly, Figure7bshows the 2D cartilage segmentation results.

Table 1 .
Eigenvalues of the Hessian matrix for different morphological structures in 3D Images.(H and L denote the high and low absolute values, respectively, and − indicates the sign of the eigenvalues.)

Table 2 .
The mean and standard deviation of Dice similarity coefficient (DSC) for different surface models.

Table 3 .
The number of iterations and running time comparisons of different surface models.

Table 4 .
The mean and standard deviation of Dice similarity coefficient (DSC) for different automatic knee cartilage segmentation methods.