Comparison of Different Approaches for Measuring Tibial Cartilage Thickness

1 Introduction

Many people experiencing knee joint pain and dysfunction suffer from osteoarthritis (OA), one of the leading causes of disability in the older population. This degenerative joint disease leads to a progressive and irreversible loss of articular cartilage [1]. To visualize the joints of affected patients, imaging modalities like ultrasound (US), magnetic resonance imaging (MRI) or contrast agent-enhanced computed tomography (CT) can be used [2], [3].

Since the mechanical properties of the knee joint behave differently under load, imaging under weight-bearing conditions is of increasing interest in order to better understand OA of articular knee cartilage [4]. This can be achieved by scanning patients in an upright standing position using a cone-beam C-arm CT (CBCT) system. A contrast agent-enhanced acquisition allows for an assessment of cartilage change in the resulting 3D reconstructions. An acquisition protocol consisting of several ten-second scans in a time interval of 15 min makes it possible to capture the progression of cartilage deformation over time.

While this 3D visualization of the joint allows for a qualitative analysis of cartilage change, a precise strain analysis requires a quantification of the articular cartilage thickness. For this purpose, a definition of a distance measure that can be applied to this complex structure is needed.

In this paper, four methods for cartilage thickness assessment are implemented and applied to segmentations of tibial bone and cartilage. The four methods include three approaches commonly used for cartilage thickness measurement (nearest neighbor, surface normals, local thickness) and an approach that has been applied to the knee joint, but not yet for measuring tibial cartilage thickness (potential field lines). The methods are chosen based on a detailed literature review (see Section Section 2). To allow for a direct comparison of the results, they are applied to the same segmentation data set, once noise-free and once deliberately corrupted by Gaussian noise. The resulting thickness maps are compared qualitatively and quantitatively to determine their differences, robustness and applicability for cartilage strain analysis.

2 Related Work

In recent years, many different approaches for measuring articular cartilage have been developed and applied. Distance measures can be performed on the original gray scale data of US, CT or MRI, but more frequently a prior segmentation of tibial and femoral bone and cartilage surfaces is required as basis for thickness computation.

A time-consuming and rather subjective approach to measuring cartilage thickness is an expert labelling in gray scale images. In studies by Kladny et al. [5], Wyler et al. [6] and Omoumi et al. [7], two raters manually measured articular cartilage thickness at well-defined positions in 2D MRI or contrast medium-enhanced CT slice images. This led to poor inter-observer agreement with only substantial correlation coefficients of 0.61–0.80 in most cases [7], relative differences of 15 % [5] and absolute differences of up to 2.51 mm [6].

In order to achieve more consistent measures, an automatic approach is beneficial. A straightforward and simple possibility to define cartilage thickness is the closest point-to-point distance between the bone and cartilage surfaces, also called nearest neighbor approach [8], [9], [10], [11]. This approach can be implemented efficiently using for example the sphere growing algorithm suggested by Li et al. [9].

Another approach to measuring cartilage thickness is to follow along the normal vector of the tibial or cartilage surface [12], [13]. However, the means of computing the vector perpendicular to a surface or point cloud vary in different studies. While Lösch et al. [14] simply defined the surface normal along the 2D gray value gradient, a more robust approach was used in recent studies by Stammberger et al. [15], [16], [17], [18], [19], namely the Euclidean Distance Transform (EDT). Other approaches utilize a curve or plane fitting to the data [20], [21] or an offset-based normal vector computation [22], [23]. Irrespective of the method for computing the normal, the cartilage thickness is subsequently defined as the Euclidean length of the vector along the normal between bone and cartilage surface.

In 1997, Hildebrand and Rüegsegger [24] defined a method for the assessment of thickness in three-dimensional data. The so-called local thickness of a volume at a point is the diameter of the largest sphere, which contains the point and which lies completely within the structure. Regarding cartilage thickness measurement, this method has initially been applied to small animal studies [25], [26]. An application for human cartilage morphology assessment can be found in a recent study by Delecourt et al. [27].

A more recent approach makes use of Laplace’s equation to define a distance measure between two surfaces [28], [29], [30], [31], [32]. The method assigns different fixed potentials to the surfaces, and by solving Laplace’s equation in the volume between, a potential field is obtained. This field can be compared to the electrostatic field of a plate capacitor [32]. Subsequently, the lengths of trajectories passing through the field define thickness. Until now, this method has been applied to measuring brain cortical structures [28], [30], patellar and femoral cartilage [29], [31] and knee joint space width [32].

Other less common approaches are a medial axis extraction for normal generation [33], midpoint sphere fitting for femoral cartilage analysis [34] or the usage of MR shape models [35], [36].

3 Materials and Methods

3.1 Data Acquisition and Labelling

The cartilage thickness measurements are applied to manual segmentations of contrast agent-enhanced CBCT reconstructions. The protocol for data acquisition was IRB-approved. The knees of one healthy subject (male, 52 years, 175 cm, 66 kg) were scanned three times using a clinical C-Arm CT system (Artis Zeego, Siemens Healthcare GmbH, Erlangen, Germany) with a flat panel detector. First, a supine CBCT acquisition was performed with the following parameters: The scan time was 20 s acquiring 496 images. The isotropic pixel resolution was 0.308 mm/pixel and the detector had a size of 1240 × 960 pixels. A dose of 1.20 μGy/frame with a peak voltage of 81 kVp was applied. After the supine scan, the subject was scanned two times in an upright weight-bearing standing position with the C-arm rotating on a horizontal trajectory [37]. One scan was acquired directly after standing up and the other one after 15 min of putting load on the joint. The acquisition parameters were the same as for the supine scan, except for a shorter scan time of 10 s acquiring 248 projections. The 3D images were reconstructed using motion compensated filtered backprojection as explained in [38] based on [39], [40]. The reconstruction results in a volume of size 512³ with an isotropic voxel spacing of 0.2 mm. Since contrast agent was only injected into the right knee, the reconstructed volumes only contained this knee. A sagittal slice of the reconstructed gray-scale image is shown in Figure 1A.

Figure 1:

Segmentation in the sagittal view of the gray-scale data. (A) Original gray scale image. (B) Overlay of an exemplary segmentation: tibia shaded in red; cartilage contact region segmented in standing scans marked as solid red line; larger region segmented in supine scan marked as dotted red line. (C) Segmentation only; region on tibial surface used for thickness calculation marked as solid black line.

Five experienced raters independently segmented the whole tibial bone and the lateral and medial tibial cartilage surface from the three reconstructed gray scale volumes using AMIRA software (AMIRA, Mercury Computer Systems, Berlin, Germany). This resulted in dense point clouds. Figure 1B shows an overlay of an exemplary segmentation in the sagittal plane. In prospect to a future analysis of cartilage change over time, the segmentations of the three scans were registered to the same coordinate system before further processing. Since the tibial bone has the same shape in all three scans, its segmentation was used to rigidly register the two standing scans to the supine scan. Afterwards, the implemented methods for thickness computation were applied to all three scans separately.

Additionally, the methods’ robustness to noise present in the data or the segmentations was evaluated. For this purpose, Gaussian noise with zero mean and unit variance was added to the z-axis of the segmented voxels of all individual segmentations.

3.2 Compared Methods

The four approaches implemented in this work are the nearest neighbor (NN) [9], surface normals (SN) [14], local thickness (LT) [24] and potential field line (FL) [32] methods. The thickness values computed by the implemented methods were always assigned to the registered tibial surfaces. This distinction is particularly important for the nearest neighbor and surface normal approach, since their results heavily depend on the starting surface of the algorithm. In contrast, local thickness and field lines work on the volume between tibia and cartilage surface. An overview is given in Figure 2.

Figure 2:

Schematic description of the four implemented methods, sagittal view. The lower bold black curved line depicts the surface of the tibial bone, the upper one the cartilage surface. Red elements explain the respective method. (A) Nearest neighbor. (B) Surface normals. (C) Local thickness. (D) Field lines.

To decrease computational effort, the segmented volumes of tibial bone were reduced to the part relevant for cartilage thickness computation. For this purpose, the region of the tibial bone situated below the cartilage was extracted. This was done by projecting the points belonging to cartilage orthogonally onto the tibia voxels along the z-axis of the volume, which approximately corresponds to the craniocaudal direction. The resulting surface was then expanded to the sides by 5 voxels respectively 1 mm to ensure that no relevant data was discarded. In the following sections, this smaller surface is referred to as “tibial surface”. This region is depicted as black solid line in the sagittal view of an example segmentation shown in Figure 1C.

3.2.1 Nearest Neighbor (NN)

Since the evaluation was not time critical and the input point clouds in most cases contained less than 10,000 points, a brute force search algorithm was implemented for the nearest neighbor method. For each point on the tibial surface, the Euclidean distance to each point on the cartilage surface is calculated. The cartilage thickness assigned to the tibial surface point is the minimum of those distances.

3.2.2 Surface Normals (SN)

As described in Section Section 2, there exist multiple approaches to computing normal vectors on surfaces. In this study, a normal computation based only on point information as described by Hoppe et al. [21] is applied. For each point on the tibial surface, a plane is fitted to its surrounding points (3 voxels to each side), and the normal vector is computed as the eigenvector belonging to the largest eigenvalue of their covariance matrix. Robustness of the normal vector computation can be increased by enlarging the number of surrounding points, while at the same time smoothing out surface information. The intersection of a line along the normal vector with the cartilage surface is computed. The starting point gets assigned the Euclidean length of this line as thickness value.

3.2.3 Local Thickness (LT)

To compute the volumetric local thickness between tibial surface and cartilage surface, an implementation by Dougherty et al. [41] developed as plugin for ImageJ [42] is used. Each voxel in the volume is assigned the diameter of the largest sphere that can be fit in the volume and that contains the point. In this way, also the voxels on the tibial surface get assigned a value that is stored as the cartilage thickness.

3.2.4 Field Lines (FL)

A more recent approach defines the thickness between two surfaces as the length of field lines passing through a potential field. This field is computed by assigning different potentials to the two surfaces and solving the Laplacian equation in the volume between. One advantage of this method is that field lines are independent of the starting surface. Furthermore, their curved path through the volume might be a better representation of the complex cartilage morphology. The computation of the field between the two surfaces is based on an iterative update described in [32]. The voxels of the two surfaces are assigned initial potentials of Φ₀ and Φ₁. The field potential between the surfaces is then computed by solving the Laplace equation:

(1)∇2Φ=0

Using the Jacobi method to solve the equation and discretizing its result, the iterative update formula applied to the volume is

(2)Φ(x,y,z)n+1=Φ(x+1,y,z)n+Φ(x−1,y,z)n+Φ(x,y+1,z)n+Φ(x,y−1,z)n+Φ(x,y,z+1)n+Φ(x,y,z−1)n6,

with x, y and z being the 3D voxel coordinates of a point, n being the current iteration, and Φ(x,y,z)n being the potential of a voxel at coordinate (x, y, z) at iteration n. The algorithm converges if the root mean squared difference of the field between two iterations is below a defined threshold of 10⁻⁹. The field lines are then generated by following the 3D gradient of the field through the volume with a small step width of 0.1 mm, until the other surface is reached. The length of the field lines is assigned as cartilage thickness.

4 Results

Figure 3 and Figure 4 show the results for all four methods of the lateral cartilage of the first standing scan based on the segmentations of one rater. While Figure 3 depicts for each method all voxels that this method assigned a value to, in Figure 4 only the voxels all methods assigned a value to are compared.

Figure 3:

Thickness values (mm) of lateral cartilage projected along the z-axis. All voxels that have been assigned a thickness value by the respective method are shown. (A) Nearest neighbor. (B) Surface normals. (C) Local thickness. (D) Field lines.

Figure 4:

Thickness values (mm) of lateral cartilage projected along the z-axis. Only the voxels that all methods assigned a value to are shown. (A) Nearest neighbor. (B) Surface normals. (C) Local thickness. (D) Field lines.

Table 1 shows the mean CC and RMSE values of the four methods over all segmentations. The mean CC is above 0.9 between all methods, and the thickness values vary on average less than 0.5 mm.

In Table 2, the mean thickness values of the four methods are shown. Presented are the mean values over all raters for each scan, split up in lateral and medial cartilage. The lateral cartilage is consistently assigned higher thickness values than the medial cartilage. For each scan and segmentation, the NN method assigned on average the lowest values, while the mean thickness was highest for LT. SN and FL lie between those two methods and achieve similar results.

Figure 5 depicts the thickness distributions already shown in Figure 3 compared to the thickness distributions resulting from noisy input data. Table 3 and Table 4 show the same metrics as Table 1 and Table 2 resulting from noisy input data. The correlation between the LT and all other methods decreases considerably and the RMSE is around 1 mm. Between the remaining methods, a slight decrease of correlations and approximately a duplication of the RMSE can be observed. The mean thickness values change only slightly for SN, are around 0.3 mm lower for NN and around 0.15 mm lower for FL and decrease markedly for the LT method.

Figure 5:

Thickness values (mm) of lateral cartilage projected along the z-axis. Top row: original segmentations as input data, bottom row: noisy segmentations as input data. (A) Nearest neighbor. (B) Surface normals. (C) Local thickness. (D) Field lines.

5 Discussion

Qualitative as well as quantitative analysis show that all four methods result in similar thickness distributions. It is noticeable that the RMSE between SN and FL approach lies on average below the voxel spacing of 0.2 mm, implying utmost similarity in this image resolution. This finding is confirmed by their results’ visual resemblance. The highest RMSE is obtained for the comparison of NN and LT, which is also reflected in the mean thickness values depicted in Table 2. Since the NN approach searches for the closest voxel on the opposite surface, this method forms a lower bound of thickness estimation. In contrast, the LT algorithm tries to fit the largest possible spheres into the volume, resulting in an upper bound for thickness calculation.

Since bone and cartilage surface in the healthy knee joint are known to be smooth [43], also the thickness calculation should result in smooth distributions. Comparing the thickness maps displayed in Figure 3 and Figure 4, this is only given for NN, LT and to some extent for FL. The SN method’s thickness distribution looks considerably less smooth than for the other methods. Furthermore, some outliers appear in the SN thickness map, see Figure 3B. These findings can be explained by the fact that this method is more sensitive to surface irregularities present in the segmentations.

The main differences between the thickness distributions can be seen at the border of the segmented cartilage, as depicted in Figure 3. Especially the consistently larger values for border pixels are striking for the FL method (Figure 3D). They arise as a result of the behavior of a potential field at the border of its source potentials, comparable to the fields of a magnet or plate capacitor. Since each rater individually defined the extent of the contact area between femoral and tibial cartilage, each segmentation of the same data set had a different size. It is therefore questionable how the methods’ differences at border regions will influence further processing of the thickness maps. Comparing only the voxels all methods assigned a value to (Figure 4), the thickness distributions over the segmented cartilages look similar, albeit varying in value.

Applying the methods to noisy data showed their robustness to possible data or segmentation inaccuracies. What can be observed from Figure 5 is that all distributions become less smooth. The method most affected by noise is LT, since adding noise generates small dips on the surface of the tibia that only spheres with small diameter fit in. The resulting distribution of the SN’s method looks similar, but a closer look shows there are plenty of outliers not only in the border regions. This confirms the method’s sensitivity to surface irregularities mentioned before. The NN and FL method are least affected by noise and result in slightly lower thickness values.

Although it is expected that the cartilage is thickest for the lying scan and decreases under load, the mean thickness values for lying scan 1 are consistently lower than for standing scan 2, see Table 2. The reason is that the raters were told to segment tibial cartilage only in the contact area between tibial and femoral cartilage for the standing scans, but to segment a larger area reaching out to thinner border regions of cartilage for the supine scan, since contact area is harder to define there. In Figure 1B and C, the region segmented for the standing scans is depicted as solid red line, whereas the larger region segmented in the supine scan is marked as red dotted line. Including the border regions to the thickness computation leads to a decrease of the mean thickness value and simultaneously to an increase of its standard deviation. For this reason, it is unsuitable for a fair comparison between supine and standing configurations. Assuming that the segmented regions for the two standing scans have approximately the same size, a tendency can be observed based on Table 2: The lateral cartilage thickness, if at all, only marginally decreases under load, while the change of medial cartilage thickness is higher with on average 0.2 mm. However, to allow for a valid comparison of the cartilage thickness distributions over time, only the overlapping regions that have been segmented in all three scans should be considered. Future work on the analysis of cartilage strain will concentrate on this issue. Since the focus of this paper was the comparison of the performance of different thickness measurement approaches, it was decided to use the complete segmentations.

6 Conclusion

A comparison of four different methods for tibial cartilage thickness computation has been presented. This comparison was conducted to examine their differences with regard to a subsequent cartilage strain analysis. The compared methods were chosen based on a comprehensive literature review. The nearest neighbor, surface normals and local thickness methods have already been applied for the purpose of measuring tibial cartilage thickness, in contrast to the field line approach that was previously used to measure femoral cartilage thickness and knee joint space width.

From the present evaluation, it can not clearly be decided which method is best suited for knee cartilage thickness assessment. Since the NN and LT methods form lower and upper bounds for the thickness, the true thickness might lie in between those values. The LT method performed markedly worse on noisy data. The SN method, that was already sensitive to variations in the data, produced an increased number of outliers under the presence of noise. Based on the limited data examined in this study, the NN and FL methods seem to be superior to the other methods. Since the FL method is, in contrast to the NN method, symmetric to the two surfaces, it might be better suited for cartilage thickness computation. However, the decision for FL and against the other methods is not definite. Further analysis on genuine segmentations of more data sets will be part of future studies and will clarify which method is most suited for the analysis of cartilage thickness.

As previously stated, the prospective goal is a computation of cartilage strain based on thickness change. Future work will concentrate on this issue and will show whether and to which degree the choice of thickness computation method influences the strain calculation.