Objective. To develop a precise three-dimensional (3D) segmentation technique for bone erosions in high-resolution peripheral quantitative CT (HR-pQCT) datasets to measure their volume, surface area and shape parameters. Assessment of bone erosions in patients with RA is important for diagnosis and evaluation of treatment efficacy. HR-pQCT allows quantifying periarticular bone loss in arthritis.
Methods. HR-pQCT scans with a spatial resolution of about 120 µm of the second to fourth metacarpophalangeal joints were acquired in patients with RA. Erosions were identified by placing a seed point in each of them. After applying 3D segmentation, the volume, surface area and sphericity of erosions were calculated. Results were compared with an approximation method using manual measurements. Intra- and interoperator precision analysis was performed for both the 3D segmentation and the manual technique.
Results. Forty-three erosions were assessed in 18 datasets. Intra- and interoperator precisions (RMSCV/RMSSD) for erosion volume were 5.66%/0.49 mm3 and 7.76%/0.76 mm3, respectively. The correlation between manual measurements and their simulation using segmentation volumes was r = 0.87. Precision errors for the manual method were 15.39% and 0.36 mm3, respectively.
Conclusion. We developed a new precise 3D segmentation technique for quantification of bone erosions in HR-pQCT datasets that correlates to the volume, shape and surface area of the erosion. The technique allows fast and effective measurement of the erosion size and could therefore be helpful for rapid and quantitative assessment of erosion size.
RA is a chronic inflammatory autoimmune disease that affects synovial joints and leads to the destruction of periarticular bone. Bone erosions are localized lesions with a break in the cortical shell. Since bone erosions are closely related to disease activity, they are an early prognostic indicator and an important clinical parameter for monitoring treatment efficacy [1, 2]. It is therefore desirable to detect them as early as possible with high precision in order to quantify small changes.
Currently, conventional radiographs (CRs) are used to semi-quantitatively assess bone erosions in patients with RA or PsA. However, due to their projectional character, the use of CRs results in an underestimation of the number and size of erosions and therefore, probably, disease activity. Furthermore, it has been reported that it takes 6–12 months before structural changes become evident on CRs , impeding the early validation of treatment efficacy. Some methods for the quantification of erosions have also been developed for CT and MRI datasets. Duryea et al. , using CT data, described a semi-automated outlining of the periosteal surface followed by an erosion segmentation based on region growing. The studies performed by Døhn et al. [5–7], also using CT data, and by Bird et al.  and Crowley et al. , using MRI data, relied on manual outlining of the erosions slice by slice. In contrast, Emond et al.  employed a true three-dimensional (3D) segmentation of erosions in MRI data. Only the placement of a seed point and the selection of one of five pre-determined parameter sets were required.
While a trained operator may produce reliable results, manual outlining can be very time consuming. Moreover, a slice-wise approach does not take the true 3D erosion structure into account. The strength of MRI is the depiction of pathological changes in soft tissue; however, the limited spatial resolution hides morphological details. While the spatial resolution of standard clinical whole-body CT scanners is better than that of MRI, it is still inadequate to depict small erosions. Recently we have applied the new technique of high-resolution peripheral quantitative computed tomography (HR-pQCT), which was originally developed for the quantification of trabecular bone architecture in the distal tibia and radius, to estimate the size of bone erosions in the MCP joints . We used simple distance measures obtained manually from selected two-dimensional (2D) slices of the HR-pQCT dataset to approximate the erosion volume [12, 13].
The aim of this study was to develop a more automated 3D erosion segmentation to improve the accuracy of size and shape measurements. In addition, we wanted to measure BMD in the erosion and its immediate vicinity.
Materials and methods
HR-pQCT datasets of 18 RA patients (16 females, 2 males, mean age 52.6 ± 15 years) from the Rheumatology Outpatient Clinic of the University of Erlangen were included in the study. All patients fulfilled the 2010 ACR/European League Against Rheumatism (EULAR) classification criteria for RA . The study was performed in accordance with the Declaration of Helsinki. Approval from the local ethics committee (Ethics Committee of the Friedrich-Alexander-University of Erlangen-Nuremberg) and the national radiation safety agency (Bundesamt für Strahlenschutz) and informed consent was obtained for the study.
High-resolution peripheral quantitative CT
All patients received an HR-pQCT scan of the MCP joints (2–4) of the dominant hand (XtremeCT, Scanco Medical AG, Brüttisellen, Switzerland) with an isotropic voxel size of (0.082 mm)3. A total number of 322 slices was acquired per dataset, 80 distally and 242 proximally relative to the top of metacarpal head 3. It is known from previous studies (e.g. Buckland-Wright ) that erosions in the MCP joints often occur close to the attachments of ligaments and at periarticular margins, so for anatomical reasons, a larger part of the scan covered the metacarpal bone.The hand was positioned in stretched posture and padded . Scan time was about 8 min.
Automated erosion quantification (AUTO3D)
The method consists of three major steps: segmentation of the periosteal bone surface, segmentation of the erosions and calculation of the quantitative parameters. All steps are integrated into the Medical Image Analysis Framework (MIAF) developed at the Institute of Medical Physics (IMP), Erlangen, Germany.
Segmentation of the periosteal bone surface
For the segmentation of the periosteal surface we applied a technique originally developed for the proximal femur . As an initial step, the metacarpal bones have to be separated from the phalangeal bones, because each bone must be segmented separately. For this purpose a sphere is placed around each metacarpal head to mask it during segmentation of the corresponding proximal phalangeal bone (Fig. 1A). The periosteal segmentation is based on volume growing with local adaptive thresholds and morphological operations. It is fully automated, but the user has the ability to correct the results.
Segmentation of erosions
The detection of erosions has not yet been automated. Erosions still have to be identified by the operator, who has to place a seed point in each of them. Then a two-stage segmentation approach is used.
First step: The level set method, a widely used segmentation technique , was used. A closed 3D surface, such as a sphere, is iteratively inflated until it stops at voxels with high grey values, resulting in a specific volume of interest (VOI), termed VOI-A. The level set surface is confined by the periosteal surface, but it may leak out into the bone marrow space if the erosion is not completely bounded by regular trabecular or newly formed sclerotic bone (Fig. 1B).
Second step: This step is only necessary once VOI-A has leaked out. It is based on the morphological concepts of ultimate erosion and skeletonization by influence zones . VOI-A is iteratively eroded with a small morphological structuring element. As the volume decreases, it will eventually separate into several pieces. These pieces are detected and recorded. This process continues until VOI-A has been completely eroded. Afterwards, each recorded individual piece is dilated separately, until all pieces join again (Fig. 2). This allows the removal of the smaller pieces, resulting in VOI-B, which defines the erosion. In cases where this process eliminates the connection to the periosteal surface, it is re-created by volume growing. Examples of segmented erosions are shown in Fig. 3.
The erosion volume V is obtained by multiplying the number of voxels of VOI-B by the voxel size. For the calculation of the surface area A, VOI-B is converted to a triangle mesh and then the areas of all the triangles are summed. The sphericity Ψ is a measure of how closely a shape resembles a perfect sphere . It can be calculated using the following formula:
BMD was calculated inside the erosion and in four concentric VOIs of increasing distance obtained by morphological dilations of VOI-B. Since the XtremeCT datasets are internally calibrated to BMD, BMD is simply the mean grey value of all voxels within a given VOI.
Manual analysis (MANUAL)
The manual approach has already been described in detail . In short, it approximates the erosion volume by constructing a half-ellipsoid from manual distance measurements. The three slices showing the largest extent of the cortical break in two perpendicular planes, i.e. along either the transverse and coronal or transverse and sagittal directions (depending on the location of the erosion), are identified by the operator. In both planes the cortical break is closed by a line and the lengths of both are measured. Perpendicular to these two lines, the corresponding maximum depths are measured (Fig. 4A). The two cortical break measurements and the greatest depth are taken to be the axes of a half-ellipsoid (Fig. 4B). Its volume is used as an approximation of the erosion volume.
Simulation of manual measurements (MANUALSIM)
We also simulated the manual approach using the erosion VOI obtained from the automated analysis. For this purpose the surface area of the cortical break Acort was determined from the intersection of the periosteal VOI and the erosion VOI. Perpendicular to each point of this surface, a ray was cast to the boundary of the erosion to find the maximum depth dmax. Acort and dmax were then used to construct the half-ellipsoid whose volume was used to approximate the erosion volume found in the manual process.
Inter- and intraoperator precision errors were determined as the root mean square (RMS) of coefficients of variation (CVRMS) and the RMS of s.d. (s.d.RMS) as described in Glüer et al. . We calculated s.d.RMS, since for small erosions, small changes lead to high coefficients of variation, while s.d. gives the absolute deviation.
We performed intra- and interoperator analysis for AUTO3D. For the intraoperator precision, the 18 HR-pQCT datasets were analysed three times by the same operator. For the interoperator precision, three different operators analysed all datasets independently. For 12 of the 18 patients, results from the MANUAL method were available, measured by two operators independently. We used these 12 datasets to compare AUTO3D to MANUAL and performed intraoperator analysis.
A total of 43 erosions (2.39 ± 1.42 per dataset) were identified. The segmentations of 18 erosions were manually edited by at least one operator. The extent of manual editing ranged from the editing of the periosteal surface in the vicinity of the cortical break to erasing VOIs that the algorithm had falsely identified as part of an erosion.
For AUTO3D, intra- and interoperator precision errors for erosion volume, surface area and sphericity are shown in Table 1. Intraoperator precision error for erosions with volumes >10 mm3 (N = 8) was 3.02%/0.92 mm3 as compared with 6.11%/0.32 mm3 for smaller erosions (N = 35). Interoperator precision error for erosions with volumes >10 mm3 (N = 9) was 5.99%/1.53 mm3 as compared with 8.27%/0.35 mm3 for smaller erosions (N = 34). Intraoperator precision error for erosions segmented fully automatically (N = 26; 2.82%/0.11 mm3) was less than for manually edited erosions (N = 17; 8.31%/0.77 mm3), and interoperator precision error for erosions segmented fully automatically (N = 25; 2.64%/0.17 mm3) was less than for manually edited erosions (N = 18; 11.58%/1.16 mm3).
|Minimum||0.22 mm3||2.25 mm2||0.33|
|Maximum||80.49 mm3||274.48 mm2||0.78|
|Mean||9.31 mm3||37.60 mm2||0.56|
|SDRMS||0.49 mm3||2.76 mm2||0.01|
|SDRMS||0.76 mm3||4.25 mm2||0.03|
|Minimum||0.22 mm3||2.25 mm2||0.33|
|Maximum||80.49 mm3||274.48 mm2||0.78|
|Mean||9.31 mm3||37.60 mm2||0.56|
|SDRMS||0.49 mm3||2.76 mm2||0.01|
|SDRMS||0.76 mm3||4.25 mm2||0.03|
Intra- and interoperator precision of the 3D segmentation method are shown for volume, surface area and sphericity. CVRMS: RMS of coefficients of variation, SDRMS: RMS of standard deviation (see Statistical analysis section).
We also calculated precision errors (CVRMS) for the four concentric VOIs for BMD measurements. Intraoperator precision errors for BMD in VOIs 1–4 (where VOI 1 is closest to the erosion) were 1.65%, 2.78%, 3.28% and 2.62%, and interoperator precision errors were 2.43%, 3.06%, 1.35% and 2.26%, respectively.
Comparison of algorithms
A subset of 12 subjects was analysed using MANUAL. A total of 22 erosions were identified and analysed using all three methods. Fig. 5 shows correlations among the different algorithms. The correlation between MANUAL and MANUALSIM was r = 0.87, while the correlations between AUTO3D and MANUAL/MANUALSIM were r = 0.61 and r = 0.72, respectively. For the same set of erosions, interoperator precision errors (CVRMS and SDRMS, two operators) for MANUAL for the calculated volumes were 15.39% and 0.36 mm3, respectively.
We present an innovative new 3D segmentation technique (AUTO3D) for bone erosions in HR-pQCT datasets of patients with RA. HR-pQCT offers excellent spatial resolution at a very low radiation exposure. It allows precise in vivo characterization of bone erosions and assessment of bone structure. However, long scan times often cause motion artefacts that impair image analysis. Also the anatomical range of the individual scan is rather limited.
Several approaches for quantifying bone erosions in MR and CT images have been developed. Many use a pure manual segmentation approach [5–9]. In high-resolution datasets, where more details of the bone structure are visible, a slice-by-slice approach often does not give satisfying results, in particular, if erosions are not completely bounded by sclerotic bone. A segmentation algorithm that accounts for the 3D nature of the erosions overcomes these limitations. Emond et al.  presented such an algorithm based on level sets and region growing for MRI datasets. However, HR-pQCT has higher spatial resolution and depicts details that are not seen on MRI. Albrecht et al.  reported that erosions <10 mm3 were not always detectable on MRI, although visible on HR-pQCT. They also found a few lesions in the MRI datasets that turned out to be very small or not visible on HR-pQCT, suggesting that the strength of MRI is to depict inflammatory tissue rather than morphological details.
Intra- and interoperator precision errors of AUTO3D were low with no or minimal user interaction. Extensive user interactions are necessary in three cases: First, for erosions that are not bounded by sclerotic bone. In such a case the level set segmentation often leaks out into the trabecular bone or into the marrow space and even the automatic correction approach described above may not fully compensate these errors. Effective 3D editing tools are provided with our implementation; however, they introduce a subjective operator interaction. Second, in case of large cortical breaks (>5 mm2), the original periosteal surface may not be restored by the periosteal segmentation algorithm and again operator interactions may be required. Third, motion artefacts, which are common with the use of the XtremeCT device due to the long scan times, may prevent an AUTO3D analysis of individual erosions or even of the complete scan. No data have been reported yet on the magnitude of the problem in patients with RA. In the default application of XtremeCT in the ultra distal radius in osteoporotic patients, where scan times are about 3 min, in a recent multicentre study, in 25% of the patients the assessment of trabecular bone structure was impaired by motion artefacts . User interactions increase precision errors, in particular, if users disagree about the exact border of an erosion, which is often the case if the boundary is not clearly identifiable.
The manual analysis illustrated in Fig. 4 is easy to use and less impacted by image quality since a human reader is more capable of interpreting blurred datasets. Also, for line measurements, only a few reference points need to be identified. An exact 3D approach depends to a greater degree on unimpaired image data than a manual approach. However, the manual analysis only weakly correlates with the RAMRIS score (r = 0.514) . The approximation of the erosion volume by a half-ellipsoid tends to underestimate the erosion volume determined by the AUTO3D segmentation, as indicated by the low correlation between the two techniques. Although no specific accuracy analysis was performed in this study, it is likely that the AUTO3D results are more accurate because they are based on 3D segmentation techniques instead of an estimate from simple line measurements. This argument is supported by the results of MANUALSIM, where the cortical break area together with the depth of the erosion was used to estimate erosion volume. This approach is similar to the MANUAL analysis. Therefore the correlation between MANUALSIM and MANUAL is relatively high and shows that the segmentation alone does not have a great impact on the MANUAL results. Instead, the difference between MANUAL and AUTO3D can be explained by the advanced analysis of the segmented erosion.
Another advantage of the AUTO3D algorithm is the possibility to determine surface area and shape parameters. The clinical relevance of the sphericity parameter has not been evaluated yet, although it has been shown that erosion shape differs among erosive diseases like RA and PsA, indicating different underlying pathophysiological concepts . Thus sphericity and other shape parameters characterizing erosions may be helpful in further research of this topic.
One limitation of the study was the relatively small number of subjects. However, the main aim was to demonstrate the feasibility and potential advantages of an automated approach to measure erosion volume with high spatial resolution. This study is a further step in establishing HR-pQCT as a diagnostic and monitoring tool in rheumatology, but obviously further clinical studies are required for its validation.
High spatial resolution imaging allows for the quantification of small bone erosions in RA.
3D image processing captures the size and shape of bone erosions in RA.
Funding: This work was supported by the Bundesministerium für Bildung und Forschung (BMBF, project ANCYLOSS).
Disclosure statement: The authors have declared no conflicts of interest.