Open AccessArticle

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

Jinke Wang

Changfa Shi

^2,3,*

Yuanzhi Cheng

⁴

Xiancheng Zhou

² and

Shinichi Tamura

⁵

Department of Software Engineering, Harbin University of Science and Technology, Rongcheng 264300, China

Mobile E-business Collaborative Innovation Center of Hunan Province, Hunan University of Commerce, Changsha 410205, China

Key Laboratory of Hunan Province for New Retail Virtual Reality Technology, Hunan University of Commerce, Changsha 410205, China

⁴

School of Computer Science and Technology, Harbin Institute of Technology, Harbin 150001, China

⁵

NBL Technovator Co., Ltd., Osaka 5900522, Japan

Author to whom correspondence should be addressed.

Symmetry 2018, 10(11), 591; https://doi.org/10.3390/sym10110591

Submission received: 26 September 2018 / Revised: 30 October 2018 / Accepted: 1 November 2018 / Published: 4 November 2018

(This article belongs to the Special Issue Information Technology and Its Applications 2021)

Download

Browse Figures

Figure 1
Flowchart of the proposed scheme on cartilage segmentation. "> 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
External forces acting upon B-spline surface points and control points. "> Figure 4
An illustration of the knee MR imaging. One 2D slice from sagittal MR sequence with anatomical structures annotated. "> Figure 5
Segmentation results of articular cartilage on one knee joint with suspected osteoarthritis (OA). (a) Original image. (b) 3D visualization of the femoral cartilage. "> Figure 6
Segmentation results of articular cartilage on one hip joint with mild OA. (a) Original image, (b) Segmented cartilage. "> Figure 7
Segmentation results of articular cartilage on one hip joint with moderate OA. (a) Original image, (b) Segmented cartilage. ">

Versions Notes

Abstract

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 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.

Keywords:

cartilage segmentation; B-spline model; vector field convolution; MRI

1. 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.

Previous studies on cartilage segmentation can be roughly divided into two categories [3]. The first category mainly relies on the pure characteristic information of the image itself, e.g., thresholding [4], region growing [5], and edge detection [6,7,8].

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.

2. 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.

2.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:

\{\begin{matrix} \frac{\partial I}{\partial t} = div (g (\nabla I) \nabla I) \\ I (t = 0) = I_{0} \end{matrix}

(1)

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 (

Δ x y

) in the projective plane. Thus, the images need to be resampled to make the two values (

Δ x y

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.

2.2. Rough Segmentation of Cartilage

2.2.1. 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:

H_{0, s} = (\begin{matrix} I_{x x} I_{x y} I_{x z} \\ I_{y x} I_{y y} I_{y z} \\ I_{z x} I_{z y} I_{z z} \end{matrix}), I_{a b} = \frac{\partial I}{\partial a \partial b}

(2)

For the Hessian matrix

H_{0, s}

, the absolute values of its three eigenvalues can be sorted in an ascending order, i.e.,

λ_{1} \leq λ_{2} \leq λ_{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]:

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.
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:

λ_{1} \approx λ_{2} \approx 0 and λ_{2} ≪ λ_{3}

(3)

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:

V_{0} (s) = \{\begin{matrix} 0 if λ_{3} > 0, \\ \exp (- \frac{{R_{A}}^{2}}{2 α^{2}}) \exp (- \frac{{R_{B}}^{2}}{2 β^{2}}) (1 - \exp (- \frac{S}{2 c^{2}})) \end{matrix}

(4)

where

R_{A} = \frac{λ_{2}}{λ_{3}}

R_{B} = \frac{λ_{1}}{\sqrt{λ_{2} λ_{3}}}

are two specially designed geometric ratios as in Reference [30],

S = \sqrt{λ_{1}^{2} + λ_{2}^{2} + λ_{3}^{2}}

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_{m i n}

s_{m a x}

V_{0} = max_{s_{m i n} \leq s \leq s_{m a x}} V_{0} (s)

(5)

2.2.2. 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.

2.3. 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.

2.3.1. 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) \in 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)

v (x, y, z) = f (x, y, z) \otimes k (x, y, z)

(6)

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:

k (x, y, z) = m (x, y, z) n (x, y, z)

(7)

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:

n (x, y, z) = \{\begin{matrix} [0, 0, 0], & if (x, y, z) = (0, 0, 0), \\ \frac{x}{\sqrt{x^{2} + y^{2} + z^{2}}}, - \frac{y}{\sqrt{x^{2} + y^{2} + z^{2}}}, - \frac{z}{\sqrt{x^{2} + y^{2} + z^{2}}}], & otherwise \end{matrix}

In addition,

m (x, y, z)

is defined as a decreasing function of the distance from the kernel origin [27]:

m (x, y, z) = \frac{1}{{(\sqrt{x^{2} + y^{2} + z^{2}} + ε)}^{γ}}

(8)

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.

2.3.2. 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:

P_{i, j} (u, w) = u M_{R} B_{R_{i j}} M_{R}^{T} w

(9)

where

u = [u^{3} u^{2} u 1]

and

w = [w^{3} w^{2} w 1]

are parameter vectors,

B_{R_{i j}}

is the

4 \times 4

control point matrix with 16 control points ranging from

Q_{i - 1, j - 1}

Q_{i + 2, j + 2}

, and

M_{R} = [\begin{matrix} - \frac{1}{6} & \frac{1}{2} & - \frac{1}{2} & \frac{1}{6} \\ \frac{1}{2} & - 1 & \frac{1}{2} & 0 \\ - \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ \frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 \end{matrix}]

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:

P_{i j} = U M_{R} B_{R_{i j}} M_{R}^{T} W^{T}

(10)

where

U = [\begin{matrix} u_{1}^{3} & u_{1}^{2} & u_{1}^{1} & 1 \\ u_{2}^{3} & u_{2}^{2} & u_{2}^{1} & 1 \\ \dots & \dots & \dots & \dots \\ u_{n}^{3} & u_{n}^{2} & u_{n}^{1} & 1 \end{matrix}], W = [\begin{matrix} w_{1}^{3} & w_{1}^{2} & w_{1}^{1} & 1 \\ w_{2}^{3} & w_{2}^{2} & w_{2}^{1} & 1 \\ \dots & \dots & \dots & \dots \\ w_{n}^{3} & w_{n}^{2} & w_{n}^{1} & 1 \end{matrix}]

, 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 \times n

surface patch matrix with surface points ranging from

P_{i j} (u_{1}, w_{1})

P_{i j} (u_{n}, w_{n})

2.3.3. 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_{i j}}

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).

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 \times n

surface points, the corresponding control points

B_{R_{i j}}

can be calculated as follows:

B_{R_{i j}} = M_{R}^{- 1} {(U^{T} U)}^{- 1} U^{T} P_{i j} W {(W^{T} W)}^{- 1} {(M_{R}^{T})}^{- 1}

(11)

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_{i j} (u, w)

of the surface is

F P_{i j} (u, w)

, then the external force that exerts on the

(i, j)

-th surface patch can be represented by an

n \times n

matrix

{F P}_{i, j}

, with elements ranging from

F P_{i j} (u_{1}, w_{1})

F P_{i j} (u_{n}, w_{n})

Assuming that

Δ B_{R_{i j}}

denotes the changes of control point matrix in each deformation step, then from Equation (11),

Δ B_{R_{i j}}

can be obtained as follows:

Δ B_{R_{i j}} = η M_{R}^{- 1} {(U^{T} U)}^{- 1} U^{T} {F P}_{i, j} W {(W^{T} W)}^{- 1} {(M_{R}^{T})}^{- 1}

(12)

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:

B_{R_{i j}}^{'} = B_{R_{i j}} + Δ B_{R_{i j}}

(13)

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.

3. Experimental Results and Analysis

3.1. 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.

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_{m i n} = 0.5

s_{m a x} = 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.

3.2. 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

D S C = \frac{2 A \cap 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.

3.3. Comparison of Different Methods

Figure 5, Figure 6 and Figure 7 show three typical cartilage segmentation results for the proposed method. Figure 5 shows 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 Figure 5a,b, respectively.

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.

4. 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.

Author Contributions

Conceptualization, C.S.; Methodology, J.W.; Software, J.W. and C.S.; Validation, C.S., Y.C., and S.T.; Formal Analysis, J.W.; Investigation, J.W. and C.S.; Resources, Y.C.; Data Curation, J.W.; Writing-Original Draft Preparation, J.W.; Writing-Review & Editing, C.S. and X.Z.; Visualization, J.W.; Supervision, C.S. and S.T.; Project Administration, S.T.; Funding Acquisition, J.W., C.S., Y.C. and X.Z.

Funding

This work was supported in part by the National Natural Science Foundation of China (No. 61741106,61701178,61571158,61471170), the National Natural Science Foundation of Hunan Province (No. 2018JJ3256), the Innovative Personnel Project of Heilongjiang Province (No. UNPYSCT-2017089), the Innovative Personnel Project of Harbin (No. RC2017QN010006), and China Scholarship Council.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

VFC	vector field convolution
OA	osteoarthritis
GVF	gradient vector flow
BCI	bone-cartilage-interface

References

Aprovitola, A.; Gallo, L. Knee bone segmentation from MRI: A classification and literature review. Biocybern. Biomed. Eng. 2016, 36, 437–449. [Google Scholar] [CrossRef]
Shan, L.; Zach, C.; Charles, C.; Niethammer, M. Automatic atlas-based three-label cartilage segmentation from MR knee images. Med. Image Anal. 2014, 18, 1233–1246. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Pham, D.L.; Xu, C.; Prince, J.L. Current Methods in Medical Image Segmentation. Ann. Rev. Biomed. Eng. 2000, 2, 315–337. [Google Scholar] [CrossRef] [PubMed]
Hardy, P.A.; Nammalwar, P.; Kuo, S. Measuring the thickness of articular cartilage from MR images. J. Magn. Resonan. Imaging 2001, 13, 120–126. [Google Scholar] [CrossRef] [3.0.CO;2-7/pdf" target='_blank'>Green Version]
Pakin, S.K.; Tamez-Pena, J.G.; Totterman, S.; Parker, K.J. Segmentation, Surface Extraction, and Thickness Computation Of Articular Cartilage; Medical Imaging 2002: Image Processing; International Society for Optics and Photonics: Bellingham, WA, USA, 2002; Volume 4684, pp. 155–167. [Google Scholar]
Mahmood, N.; Shah, A.; Waqas, A.; Abubakar, A.; Kamran, S.; Zaidi, S.B. Image segmentation methods and edge detection: An application to knee joint articular cartilage edge detection. J. Theor. Appl. Inf. Technol. 2015, 71, 87–96. [Google Scholar]
Jaremko, J.; Cheng, R.; Lambert, R.; Habib, A.; Ronsky, J. Reliability of an efficient MRI-based method for estimation of knee cartilage volume using surface registration. Osteoarthr. Cartil. 2006, 14, 914–922. [Google Scholar] [CrossRef] [PubMed]
Carballido-Gamio, J.; Bauer, J.S.; Stahl, R.; Lee, K.Y.; Krause, S.; Link, T.M.; Majumdar, S. Inter-subject comparison of MRI knee cartilage thickness. Med. Image Anal. 2008, 12, 120–135. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Ghosh, S.; Beuf, O.; Ries, M.; Lane, N.E.; Steinbach, L.S.; Link, T.M.; Majumdar, S. Watershed segmentation of high resolution magnetic resonance images of articular cartilage of the knee. In Proceedings of the 22nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Cat. No. 00CH37143), Chicago, IL, USA, 23–28 July 2000; Volume 4, pp. 3174–3176. [Google Scholar]
Grau, V.; Mewes, A.U.J.; Alcani, M.; Kikinis, R.; Warfield, S.K. Improved watershed transform for medical image segmentation using prior information. IEEE Trans. Med. Imaging 2004, 23, 447–458. [Google Scholar] [CrossRef] [PubMed]
Folkesson, J.; Dam, E.B.; Olsen, O.F.; Pettersen, P.C.; Christiansen, C. Segmenting articular cartilage automatically using a voxel classification approach. IEEE Trans. Med. Imaging 2007, 26, 106–115. [Google Scholar] [CrossRef] [PubMed]
Heimann, T.; Meinzer, H.P. Statistical shape models for 3D medical image segmentation: A review. Med. Image Anal. 2009, 13, 543–563. [Google Scholar] [CrossRef] [PubMed]
Iglesias, J.E.; Sabuncu, M.R. Multi-atlas segmentation of biomedical images: A survey. Med. Image Anal. 2015, 24, 205–219. [Google Scholar] [CrossRef] [PubMed] [Green Version]
McInerney, T.; Terzopoulos, D. Deformable models in Med. Image Anal.: A survey. Med. Image Anal. 1996, 1, 91–108. [Google Scholar] [CrossRef]
Lee, J.G.; Gumus, S.; Moon, C.H.; Kwoh, C.K.; Bae, K.T. Fully automated segmentation of cartilage from the MR images of knee using a multi-atlas and local structural analysis method. Med. Phys. 2014, 41, 092303. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Glocker, B.; Komodakis, N.; Paragios, N.; Glaser, C.; Tziritas, G.; Navab, N. Primal/dual linear programming and statistical atlases for cartilage segmentation. In Proceedings of the Medical Image Computing and Computer-Assisted Interventation (MICCAI’2007), Brisbane, Australia, 29 October–2 November 2007; pp. 536–543. [Google Scholar]
Solloway, S.; Hutchinson, C.E.; Waterton, J.C.; Taylor, C.J. The use of active shape models for making thickness measurements of articular cartilage from MR images. Magn. Reson. Med. 1997, 37, 943–952. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Cootes, T.F.; Taylor, C.J.; Cooper, D.H.; Graham, J. Active shape models-their training and application. Comput. Vis. Image Underst. 1995, 61, 38–59. [Google Scholar] [CrossRef]
Ye, T.; Cui, X.; Kim, H. Fully automated segmentation of cartilage from magnetic resonance images using improved 3D shape context and active shape model. Osteoarthr. Cartil. 2015, 23, A301–A302. [Google Scholar] [CrossRef]
Yin, Y.; Zhang, X.; Williams, R.; Wu, X.; Anderson, D.D.; Sonka, M. LOGISMOS—layered optimal graph image segmentation of multiple objects and surfaces: Cartilage segmentation in the knee joint. IEEE Trans. Med. Imaging 2010, 29, 2023–2037. [Google Scholar] [CrossRef] [PubMed]
Lee, S.; Park, S.H.; Shim, H.; Yun, I.D.; Lee, S.U. Optimization of local shape and appearance probabilities for segmentation of knee cartilage in 3-D MR images. Comput. Vis. Image Underst. 2011, 115, 1710–1720. [Google Scholar] [CrossRef]
Montagnat, J.; Delingette, H. A review of deformable surfaces: Topology, geometry and deformation. Image Vis. Comput. 2001, 19, 1023–1040. [Google Scholar] [CrossRef]
Stammberger, T.; Eckstein, F.; Michaelis, M.; Englmeier, K.H.; Reiser, M. Interobserver reproducibility of quantitative cartilage measurements: Comparison of B-spline snakes and manual segmentation. Magn. Reson. Imaging 1999, 17, 1033–1042. [Google Scholar] [CrossRef]
Lynch, J.A.; Zaim, S.; Zhao, J.; Stork, A.; Peterfy, C.G.; Genant, H.K. Cartilage Segmentation of 3D MRI Scans of the Osteoarthritic Knee Combining User Knowledge and Active Contours; Medical Imaging 2000: Image Processing; International Society for Optics and Photonics: Bellingham, WA, USA, 2000; Volume 3979, pp. 925–936. [Google Scholar]
Kauffmann, C.; Gravel, P.; Godbout, B.; Gravel, A.; Beaudoin, G.; Raynauld, J.P.; Martel-Pelletier, J.; Pelletier, J.P.; de Guise, J.A. Computer-aided method for quantification of cartilage thickness and volume changes using MRI: Validation study using a synthetic model. IEEE Trans. Biomed. Eng. 2003, 50, 978–988. [Google Scholar] [CrossRef] [PubMed]
Fripp, J.; Crozier, S.; Warfield, S.K.; Ourselin, S. Automatic segmentation and quantitative analysis of the articular cartilages from magnetic resonance images of the knee. IEEE Trans. Med. Imaging 2010, 29, 55–64. [Google Scholar] [CrossRef] [PubMed]
Li, B.; Acton, S.T. Active Contour External Force Using Vector Field Convolution for Image Segmentation. IEEE Trans. Image Process. 2007, 16, 2096–2106. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Weickert, J. Anisotropic Diffusion in Image Processing; Teubner: Stuttgart, Germany, 1998; Volume 1. [Google Scholar]
Maymon, S.; Oppenheim, A.V. Sinc interpolation of nonuniform samples. IEEE Trans. Signal Process. 2011, 59, 4745–4758. [Google Scholar] [CrossRef]
Frangi, A.F.; Niessen, W.J.; Vincken, K.L.; Viergever, M.A. Multiscale vessel enhancement filtering. In Proceedings of the Medical Image Computing and Computer-Assisted Interventation (MICCAI’98), Cambridge, MA, USA, 11–13 October 1998; pp. 130–137. [Google Scholar]
Liao, P.S.; Chen, T.S.; Chung, P.C. A fast algorithm for multilevel thresholding. J. Inf. Sci. Eng. 2001, 17, 713–727. [Google Scholar]
Xu, C.; Prince, J.L. Snakes, Shapes and Gradient Vector Flow. IEEE Trans. Image Process. 1998, 7, 359–369. [Google Scholar] [PubMed]
Chen, X. 3D Object Segmentation Using Deformable Models. Ph.D. Thesis, Nanyang Technological University, Singapore, 2007. [Google Scholar]
Yushkevich, P.A.; Piven, J.; Cody Hazlett, H.; Gimpel Smith, R.; Ho, S.; Gee, J.C.; Gerig, G. User-Guided 3D Active Contour Segmentation of Anatomical Structures: Significantly Improved Efficiency and Reliability. Neuroimage 2006, 31, 1116–1128. [Google Scholar] [CrossRef] [PubMed]
Crum, W.R.; Camara, O.; Hill, D. Generalized Overlap Measures for Evaluation and Validation in Medical Image Analysis. IEEE Trans. Med. Imaging 2006, 25, 1451–1461. [Google Scholar] [CrossRef] [PubMed]
Li, C.; Huang, R.; Ding, Z.; Gatenby, C.; Metaxas, D.N.; Gore, J.C. A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogeneities With Application to MRI. IEEE Trans. Image Process. 2011, 20, 2007–2016. [Google Scholar] [PubMed] [Green Version]

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

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. External forces acting upon B-spline surface points and control points.

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

Figure 5. Segmentation results of articular cartilage on one knee joint with suspected osteoarthritis (OA). (a) Original image. (b) 3D visualization of the femoral cartilage.

Figure 6. Segmentation results of articular cartilage on one hip joint with mild OA. (a) Original image, (b) Segmented cartilage.

Figure 7. Segmentation results of articular cartilage on one hip joint with moderate OA. (a) Original image, (b) Segmented cartilage.

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.)

Structure	Eigenvalues			Typical Case
Structure	$λ_{1}$	$λ_{2}$	$λ_{3}$	Typical Case
Plate (Bright)	L	L	H−	Cartilage
Tubular (Bright)	L	H−	H−	Vessel
Spherical (Bright)	H−	H−	H−	Tumor

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

Methods	Cartilage
Methods	Tibial	Femoral
Traditional Model	0.774 ± 0.036	0.743 ± 0.052
Ours	0.827 ± 0.041	0.792 ± 0.043

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

Methods	100 Slices		150 Slices		200 Slices
Methods	Iterations	Time (s)	Iterations	Time (s)	Iterations	Time (s)
Traditional Model	28	108	60	192	90	287
Ours	20	64	46	108	66	174

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

Methods	Cartilage
Methods	Tibial	Femoral
Folkesson et al. [11]	0.812 ± 0.054	0.756 ± 0.042
Yin et al. [20]	0.820 ± 0.043	0.761 ± 0.061
Lee et al. [21]	0.824 ± 0.038	0.776 ± 0.055
Ours	0.827 ± 0.041	0.792 ± 0.043

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wang, J.; Shi, C.; Cheng, Y.; Zhou, X.; Tamura, S. Vector Field Convolution-Based B-Spline Deformation Model for 3D Segmentation of Cartilage in MRI. Symmetry 2018, 10, 591. https://doi.org/10.3390/sym10110591

AMA Style

Wang J, Shi C, Cheng Y, Zhou X, Tamura S. Vector Field Convolution-Based B-Spline Deformation Model for 3D Segmentation of Cartilage in MRI. Symmetry. 2018; 10(11):591. https://doi.org/10.3390/sym10110591

Chicago/Turabian Style

Wang, Jinke, Changfa Shi, Yuanzhi Cheng, Xiancheng Zhou, and Shinichi Tamura. 2018. "Vector Field Convolution-Based B-Spline Deformation Model for 3D Segmentation of Cartilage in MRI" Symmetry 10, no. 11: 591. https://doi.org/10.3390/sym10110591

APA Style

Wang, J., Shi, C., Cheng, Y., Zhou, X., & Tamura, S. (2018). Vector Field Convolution-Based B-Spline Deformation Model for 3D Segmentation of Cartilage in MRI. Symmetry, 10(11), 591. https://doi.org/10.3390/sym10110591

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Methods

2.1. Image Preprocessing

2.2. Rough Segmentation of Cartilage

2.2.1. Cartilage Enhancement

2.2.2. Image Binarization

2.3. Cartilage Refinement Based on B-Spline Deformation Model

2.3.1. External Force Field of B-Spline Deformation Model

2.3.2. The Definition of B-Spline Surface Representation

2.3.3. Deformation Process of B-Spline Deformation Model

3. Experimental Results and Analysis

3.1. Image Data and Parameters

3.2. Comparison of Different Surface Models

3.3. Comparison of Different Methods

4. Discussion and Conclusions

Author Contributions

Funding

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI