Open-source image analysis tool for the identification and quantification of cortical interruptions and bone erosions in high-resolution peripheral quantitative computed tomography images of patients with rheumatoid arthritis

Identification of bone erosions and quantification of erosion volume is important for rheumatoid arthritis diagnosis, and can add important information to evaluate disease progression and treatment effects. High-resolution peripheral quantitative computed tomography (HR-pQCT) is well suited for this purpose, however analysis methods are not widely available. The purpose of this study was to develop an open-source software tool for the identification and quantification of bone erosions using images acquired by HR-pQCT. The collection of modules, Bone Analysis Modules (BAM) – Erosion, implements previously published erosion analysis techniques as modules in 3D Slicer, an open-source image processing and visualization tool. BAM includes a module to automatically identify cortical interruptions, from which erosions are manually selected, and a hybrid module that combines morphological and level set operations to quantify the volume of bone erosions. HR-pQCT images of the second and third metacarpophalangeal (MCP) joints were acquired in patients with RA (XtremeCT, n = 14, XtremeCTII, n = 22). The number of cortical interruptions detected by BAM-Erosion agreed strongly with the previously published cortical interruption detection algorithm for both XtremeCT (r 2 = 0.85) and XtremeCTII (r 2 = 0.87). Erosion volume assessment by BAM-Erosion agreed strongly (r 2 = 0.95) with the Medical Image Analysis Framework. BAM-Erosion provides an open-source erosion analysis tool that produces comparable re-sults to previously published algorithms, with improved options for visualization. The strength of the tool is that it implements multiple image processing algorithms for erosion analysis on a single, widely available, open-source platform that can accommodate future updates.


Introduction
Rheumatoid arthritis (RA) patients can develop irreversible structural joint damage, including bone erosions, as a result of inflammatory processes [1].Erosions are characterized by an interruption in the cortical bone (i.e., discontinuities in the cortical bone that extend from the periosteum to the marrow) accompanied by underlying trabecular bone loss [1,2], and are associated with disease severity as well as declines in physical function [3,4].In clinical practice, erosions are assessed with conventional radiographs or ultrasound.For research and clinical trials, magnetic resonance imaging (MRI) and computed tomography (CT) have improved the sensitivity to detect smaller erosions; however the most sensitive modality to detect and monitor longitudinal change of erosion is high-resolution peripheral quantitative computed tomography (HR-pQCT) [5][6][7][8][9].
Erosions are a subset of cortical interruptions.Cortical interruptions, defined as a break in cortical bone with underlying trabecular bone loss, are a frequent observation with high resolution imaging [9].They are present in greater quantities [10] and more frequently show significant changes in volume [11] in RA patients than in healthy controls.Not all cortical interruptions, are considered pathological interruptions (i.e., erosions).Some cortical interruptions represent physiological vascular channels [12] and are presumed to be linear in shape to differentiate them from erosions [13].These vascular channels are hypothesized to allow for osteoclastic-mediated joint destruction and the later appearance of erosions [14][15][16].Erosions are more specific for RA patients than cortical interruptions, but they also constitute a feature of age-related structural change in healthy participants [17].Therefore, a clinical decision is still required to determine whether an erosion is related to RA.
Recently several software tools have been developed for the identification and quantification of cortical erosions and cortical interruption.Peters et al. designed a precise algorithm to automatically detect cortical interruptions (CIDA) with underlying trabecular bone loss [18,19].The algorithm is currently implemented in the proprietary Image Processing Language (IPL), commercially available through the HR-pQCT scanner manufacturer (Scanco Medical AG).In this software environment, it is challenging to visualize interruptions in all three planes, and therefore assess which interruptions correspond to true pathologic erosions.
Another approach is pursued by MIAF (Medical Image Analysis Framework), whereby an operator places a seed point where they identify an erosion.Then an automated 3D segmentation is used to determine shape and volume of the erosion.If necessary an operator can correct the segmentation [20,21].MIAF replaces an earlier manual approach to estimate erosion volume as a half-ellipsoid based on the maximum width and depth [7,13,22,23].While CIDA and MIAF are proprietary, access to MIAF is further limited as it is neither in the public domain or commercially available.
To overcome limitations in accessibility, we developed an opensource tool, Bone Analysis Module (BAM) -Erosion.Using Python and several of its image processing libraries as a language, we translated the CIDA and MIAF algorithms and integrated them as a module in 3D Slicer [24], a customizable image computing platform.Our new module currently performs the following functions using HR-pQCT images: 1) to identify cortical interruptions and underlying trabecular bone loss based on Peters et al.'s automatic method [18,19], 2) to quantify erosion volume in erosions identified by an operator by a hybrid technique involving morphological and level set operations, similar to the MIAF method, and 3) image registration to track longitudinal changes in erosion volume.In this study we quantified the agreement of BAM -Erosion with CIDA for assessment of the number and size of cortical interruptions, and MIAF for assessment of erosion volume.

Study design
Image data were obtained from single time points of previously published longitudinal studies which were approved by the Conjoint Health Research Ethics Board at the University of Calgary (REB15-0582 and REB 13-0743) [5,25].All participants provided written informed consent prior to participation.For comparison of the cortical interruption detection algorithms, scans of the metacarpophalangeal (MCP) joints by first generation HR-pQCT (XtremeCT, Scanco Medical, Brüttisellen, Switzerland) were selected from a study of patients (n = 14) who had experienced inflammatory arthritis symptoms for less than one year [13].
To evaluate the performance of the cortical interruption detection algorithm with second generation HR-pQCT (XtremeCT II, Scanco Medical) as well as erosion volume analysis, MCP scans from 22 participants were selected from a study of RA patients who had inadequate response to previous treatments and were undergoing new biologic disease-modifying antirheumatic drug therapies [14].Of these 22, 17 with confirmed erosions were utilized for a comparison of erosion volume quantification algorithms [14].For the erosion volume analysis, erosions were identified according to the Study grouP for xtrEme Computed Tomography in RA (SPECTRA) definition: a cortical interruption that spans at least two consecutive slices, in at least two orthogonal planes with underlying trabecular bone loss, and is nonlinear in shape [13].Erosion volume quantification was conducted by the same operator (SCB) for both BAM-Erosion and MIAF.

HR-pQCT image acquisition
Briefly, in both datasets, HR-pQCT scans of the 2nd and 3rd MCP joints of the participant's dominant hand were acquired, as per the SPECTRA measurement protocol [2].Standard image acquisition parameters were used to acquire 3 image stacks: a 27.0 mm long region for XtremeCT (XCT, 59.4 kVp, 900 mA, 100 ms integration time, 82.0 μm nominal isotropic resolution) or a 30.6 mm long region on XtremeCTII (XCTII, 68 kVp, 1470 mA, 43 ms integration time, 60.7 μm nominal isotropic resolution).All scans were evaluated for motion using the manufacturer's standard scoring system from 1 to 5, scans with a stack that included a score > 3 were excluded [26].For this analysis, scans were additionally excluded where cortical interruptions were misidentified because of a stack shift artifact.

HR-pQCT image pre-processing
Image preprocessing was performed for CIDA and BAM-Erosion, in which extraction of the 2nd or 3rd MCP joint was conducted using Image Processing Language (IPL v5.42).Bone masks were generated by identifying the outer, periosteal surface of the 2nd and 3rd metacarpal and phalangeal bones using an automated method [27,28].

Automated cortical interruption detection and underlying trabecular bone loss
The image processing pipeline was the same for CIDA and BAM.In order to use the same size of the structuring elements of morphological operators, the number of voxels of the elements was adjusted for scanner generation.Briefly, IPL-generated periosteal bone masks were applied.For XCT images, the mineralized bone was extracted using the manufacturer's standard evaluation protocol, which uses a Laplace-Hamming filter (laplace epsilon 0.5, cut off frequency 0.4, threshold 400).For XCTII images, the mineralized bone was segmented following standard gaussian filter (σ = 0.8) and thresholding (686 HU).The endosteal surface was defined as the surface with a constant distance of 0.3 mm (4 voxels in XCT, and 5 voxels XCTII) from the periosteal surface.BAM identifies cortical interruptions of at least 0.4 mm in diameter (5 voxels in XCT, and 7 voxels in XCTII) that are connected to both the periosteal and endosteal surfaces.For images of different spatial resolution, BAM-Erosion allows for adjustment of the parameter settings for distance between periosteal and endosteal surfaces and minimum diameter of cortical interruptions.
To determine the underlying trabecular bone loss at each cortical interruption (i.e., the cortical interruption volume, which may include vascular channels), the voids in the bone are first labeled as separate objects.Each void is eroded by a user specified distance (0.4 mm) to remove connectivity with the inner trabecular bone.Afterwards, the voids are matched to the corresponding cortical interruptions using a connectivity filter, then restored to their original size through dilation.This procedure identifies the subset of cortical interruptions and eliminates those that are not accompanied by underlying trabecular bone loss.
The number of cortical interruptions in each joint from the first dataset was determined by both the BAM-Erosion and CIDA.We additionally compared the volume of underlying trabecular bone loss for interruptions that were similarly identified by both algorithms.For the purposes of comparison in this study, all input parameters were held constant for both methods.

Table 1
The number of cortical interruptions and interruption volume observed using CIDA and BAM-Erosion on first (XCT) and second (XCTII) generation HR-pQCT scanners.The sample size (n) indicates the number of bones included for each scanner.a For XCT, 7 bones were excluded due to motion artifact, 1 due to stack artifact.For XCTII, 6 bones were excluded due to motion artifact, 24 due to stack artifact.BMI = body mass index, MC2 = 2nd metacarpal, MC3 = 3rd metacarpal, P2 = 2nd phalange, P3 = 3rd phalange, vdHSS = van der Heijde modified Sharp score [37].
M. Zhao et al.

Semi-automated erosion volume quantification
Erosions from 17 participants (corresponding to 39 bones), previously evaluated in MIAF, were selected for comparison against BAM-Erosion.Seed points were visually matched to confirm selection of the same erosions in both analyses.The same operator (SCB) manually corrected erosions in both software packages.Analyses and manual corrections with each package were conducted independently, approximately 2 years apart.Erosion volume of these manually identified erosions was defined as the voxel count of the resulting segmentation multiplied by the voxel volume.
BAM-Erosion uses a hybrid technique that combines morphological and level set operations to segment and quantify erosions.First, the periosteal surface is segmented (Fig. 1b) in a similar manner to MIAF (de-noise, threshold followed by morphological operations) [20].This approach inherits the seed-based semi-automatic approach from MIAF; the operator places a seed point inside each erosion (Fig. 1c).
In BAM-Erosion a number of image filters are applied sequentially to segment the erosions.In particular, a Gaussian filter with a user specified sigma value de-noises the image (Fig. 1d); a threshold filter selects the cortical interruptions and the trabecular voids within the periosteal surface (Fig. 1e); a distance transformation filter selects the segmented regions that are wider than a user-specified diameter (Fig. 1f); a connectivity filter (Fig. 1h) is applied in between morphological erosion (Fig. 1g) and dilation (Fig. 1i) to remove undesired trabecular voids that are not connected to the seed points; and a level set region growing filter further expands the erosion segmentation (Fig. 1j) using the squared distance map of the erosion segmentation as the speed function.The additional morphological (remove connections between voids in trabecular bone) and connectivity filter (maintain only voids connected to the seed point) were added to prevent the voids from leaking too far into the trabecular bone, and to minimize the number of manual Where location of interruptions was not congruent between CIDA and BAM-Erosion, volumes were not compared.Thus, fewer participants and interruptions were included in the volume analysis.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)corrections.BAM-Erosion uses 3D Slicer's built-in drawing tools, allowing for manual correction of the resulting segmentation.
In MIAF, erosion segmentation applies a level set algorithm directly to the manually positioned seed point without the need to set any further parameters.If required, leakage of the erosion into the trabecular bone is prevented by morphological operations based on ultimate erosion and skeletonization by influence zones [20].A variety of sophisticated 3D editing tools is available to the operator [33].
While not quantitatively evaluated here, BAM-Erosion includes a module for automated image registration.This module is an implementation of the registration method from SimpleITK.Similar to the 3D registration used in MIAF, the image registration ensures the original cortical surface, which is destroyed by the erosion, is as similar as possible in both images in order to accurately determine the change in erosion volume over time [34].

Statistical analysis
For visual comparison, cortical interruptions and erosion masks output from each software implementation were overlaid in 3D Slicer.Data organization, analysis and visualization were completed in R (v. 4.0.2) and R Studio (v.1.3.1093).The number of cortical interruptions per joint and the volume of individual erosions were described as median (interquartile range [IQR]).Bland-Altman and correlations plots were used to describe agreement between methods.Root-mean square standard deviation (RMSSD) was calculated for comparison with previously reported precision errors [35].Significance was set at an alpha level of 0.05.

BAM-Erosion compared to CIDA: cortical interruptions
The XCT dataset used for the analysis of cortical interruptions came from patients (mean age 48 (SD 14) years, 64 % female) who had a mean disease duration of 3.1 (SD 5.0) years (Table 1).The median number of cortical interruptions per bone was 3 (IQR 2-5) for CIDA and 3 (IQR 1.75-5) for BAM-Erosion (Table 1).The number of cortical interruptions per joint was strongly correlated between the two implementations (r 2 = 0.85, p < 0.001, Fig. 2).
Bland-Altman plots demonstrate good agreement between implementations.On average BAM-Erosion detected 0.72 fewer interruptions per bone than CIDA, with a 95 % limit of agreement of 3.1 interruptions; and a RMSSD of 1.34.An example of a cortical interruption is shown in Fig. 3.The mean underlying interruption volume was 2.08 (SD 3.31) mm 3 for CIDA and 2.01 (SD 2.84) mm 3 for BAM-Erosion.Again, results were strongly correlated between the two implementations (r 2 = 0.96, p < 0.001).On average, interruption volumes measured by BAM-Erosion were 0.07 mm 3 smaller than measured by CIDA with a 95 % limit of agreement of 1.47 mm 3 and a RMSSD of 0.53 mm 3 (Fig. 2).
The XCTII dataset used for the same analysis came from 22 patients (mean age 53 (SD 13) years, 77 % female) who had a mean disease duration of 11 (SD 9) years (Table 1).The median number of cortical interruptions per bone was 4 (IQR 2 to 6) for CIDA and 3 (IQR 2 to 5) for BAM-Erosion (Table 1).The number of cortical interruptions was strongly correlated between the two implementations (r 2 = 0.87, p < 0.001, Fig. 4).On average BAM-Erosion detected 0.91 fewer interruptions per bone than CIDA, with a 95 % limit of agreement of 2.36 and RMSSD of 1.06 interruptions/bone.The mean underlying interruption volume was 1.68 (SD 2.83) mm 3 for CIDA and 1.70 SD (3.17) mm 3 for BAM-Erosion.Again, they were highly correlated (r 2 = 0.93, p < 0.001).On average, BAM-Erosion found interruption volumes 0.03 mm 3 larger than CIDA with a 95 % limit of agreement of 1.74 mm 3 .An example of the bone loss detected by CIDA and BAM-Erosion is shown in Fig. 5.
An additional finding was that cortical interruption detection was highly sensitive to the bone mask generated during pre-processing.Masks were generated with an automated algorithm in CIDA and visualized in the axial plane using the Scanco platform to confirm and correct if necessary.However, upon later visualization in the coronal and sagittal planes in 3D Slicer, we noted that many of the masks did not accurately capture bone surfaces in these planes and led to erroneous detection of cortical interruptions (Supplementary Fig. 1).

BAM-erosion compared to MIAF: erosion volume
HR-pQCT images included in the erosion volume analysis came from 17 RA patients with confirmed erosions.Thirty-nine erosions from 18 joints were included in the analysis.
The median volume of the erosions was 6.7 (IQR 2.2 to 26.1) mm 3 for MIAF and 6.5 (IQR 1.6 to 19.7) mm 3 for BAM-Erosion.The erosion volume assessed by BAM-Erosion strongly correlated with MIAF (r 2 = 0.95, p < 0.001, Fig. 6).Excluding erosions with the two highest volumes reduced the correlation coefficient to r 2 = 0.83 (p < 0.001).The mean difference between the erosion volume with MIAF and BAM-Erosion was 2.3 mm 3 , the RMSSD was 4.3 mm 3 , with a 95 % limit of agreement of 14.2 mm 3 .Qualitatively, there were differences in how each implementation defined the virtual cortical surface (Fig. 5) and extension of the erosion into the trabecular bone (Supplementary Fig. 2).

Image registration
We qualitatively evaluated the performance of the image registration module.The module was able to align follow-up images to the baseline image space.The module provides several visualization options (Supplementary Fig. 3) to confirm the alignment of the images.

Discussion
In this study we assessed the performance of a newly developed tool to detect cortical interruptions and measure the volume of erosions in HR-pQCT images.The new tool, BAM-Erosions is available as a module of 3D Slicer.Overall, there was good agreement between BAM-Erosion and existing tools for cortical interruption detection and erosion volume assessment.The cortical interruption volume accounts for discontinuities in the cortex that might be vascular channels or pathological erosions, whereas the erosion volume only accounts for a subset of cortical interruptions that were identified as pathological erosions.BAM-Erosion has the advantages of (1) combining currently available erosion analyses in one tool (2) being developed in an open-source environment (python libraries and 3D Slicer), and (3) providing excellent visualization options already available in 3D Slicer.As Bland-Altman plots illustrated no dependence of cortical interruption or erosion volume on size, the results should be applicable to any patient population with erosions.Further, the BAM-Erosion modules can be modified and new functions added, as new algorithms (e.g., fully automated erosion analysis) or additional steps (e.g., joint space width) are developed.
For cortical interruption detection, the two automated implementations BAM-Erosion and CIDA showed good agreement in the number of cortical interruptions per joint.The number of interruptions (3 vs. 3.1) and the RMSSD was similar (1.2) to a previously reported scan-rescan reproducibility study (1.5) [18].However, both implementations are highly sensitive to motion and stack artifacts.Nonetheless, future studies may need to manually discard erroneous cortical interruptions, implement algorithms to automatically detect the presence of the stack artifact or acquire data with overlapping stacks [36].This study also compared erosion volume quantification between BAM-Erosion and MIAF, through the manual placement of erosion seed points.Although only two erosions fell outside the 95 % limits of  agreement, the RMSSD (4.3 mm 3 ) was higher than previously reported intra-rater precision errors for MIAF (0.49 to 2.52 mm 3 ) [5,20].When defining the borders of an erosion, BAM-Erosion and MIAF use different concepts which, not surprisingly, results in different erosion volumes (Supplementary Fig. 2).The differences can be attributed to several sources.(1) Different segmentation approaches, in which MIAF, as opposed to BAM-Erosion, virtually restores the cortical surface (Fig. 6).
(2) A difference in the method to determine the extension of the erosion in the trabecular bone.These results highlight the lack of consensus on how to define the extension of an erosion in the trabecular bone compartment, which can be further complicated in larger erosions irregular shape.(3) Finally, while the analyses with BAM-Erosion and MIAF were performed by the same operator (SCB), they were acquired 2 years apart and manual corrections may have been inconsistent between the two analyses.Despite the differences between BAM-Erosion and MIAF, the high correlation of erosion volume indicates that difference may be relatively small on average, and illustrate important challenges in the assessment of erosions from HR-pQCT images.While we did not compare BAM-Erosion and MIAF in the XCT dataset, the sources of error would remain the same, and are not dependent on voxel size, thus we expect similar performance of BAM-Erosion on XCT scans.
This study shows important challenges in the assessment of erosions from HR-pQCT images.Despite excellent spatial resolution and the ability to determine erosion volume with much higher accuracy than other in vivo imaging modalities, currently, expert knowledge is required for all techniques to decide which erosions are pathologic.However, identification of a cortical interruption as an erosion has acceptable (81 %) but not perfect agreement between expert readers [23].Manual identification of the erosions can be more tolerant against motion artifacts if they don't affect the local environment of the erosion to be segmented.Further, vascular channels (i.e., interruptions not considered to be erosions) identified in histologic specimens demonstrate complex shapes that may not be captured by the SPECTRA definition [12] and have been implicated in osteoclastic activity and synovial hyperplasia with exposure to inflammation in animal models of arthritis suggesting that they may also be relevant to pathogenesis [14][15][16].Nonetheless, fully automated analyses remain a challenge.Finally, we did not compare erosion metrics with clinical outcomes here.The ease of this analysis tool should enable future studies with larger sample sizes to investigate the association of erosions and cortical interruptions with prognostic and patient-reported outcomes.
In summary, the newly-developed, open-source, BAM-Erosion module, and its incorporation into an accessible platform (3D Slicer) enables the rapid and effective assessment of cortical interruptions and erosion volumes.In addition to good agreements with previously established algorithms, the main advantage of BAM-Erosion is its accessibility.This module will enable the opportunity, with experts and new researchers, to collaboratively develop more specific and automated techniques for erosion analysis and other bone morphological tools.Thus, the accessibility of BAM-Erosion may lead to efforts that will not only better identify cortical interruptions that are associated with bone damage but identify those associated with erosions or bone changes that may contribute disease progression.

Fig. 1 .
Fig. 1.Representation of the steps implemented by the erosion segmentation algorithm to mimic semi-automated erosion analysis performed with MIAF.

Fig. 2 .
Fig. 2. Cortical interruption detection algorithm for XCT.(A) Regression plot comparing the number of cortical interruptions observed with the CIDA and BAM-Erosion implementations.Solid line indicates regression line, dashed blue line indicates the line of unity.(B) Bland-Altman plot comparing the differences between CIDA and BAM-Erosion implementations in the number of cortical interruptions for each joint.The mean difference is shown as a solid black line and 95 % limits of agreement are shown in dashed black lines.n = 48 bones from 14 RA participants for counts of cortical interruption.(C) Regression plot comparing cortical interruption volume between implementations.(D) Bland-Altman plot comparing in cortical interruption volume between implementations.n = 66 interruptions from 12 participants for cortical interruption volume.Where location of interruptions was not congruent between CIDA and BAM-Erosion, volumes were not compared.Thus, fewer participants and interruptions were included in the volume analysis.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Visualization of cortical interruption and underlying trabecular bone loss detected by CIDA and BAM-Erosion in an XCT image of the 2nd metacarpal.

Fig. 4 .
Fig. 4. Cortical break interruption algorithm for XCTII.(A) Regression plot comparing the number of cortical interruptions observed with the CIDA and BAM-Erosion implementations.Solid line indicates regression line, dashed blue line indicates the line of unity.(B) Bland-Altman plot comparing the differences between CIDA and BAM-Erosion implementations in the number of cortical interruptions for each joint.The mean difference is shown as a solid black line and 95 % limits of agreement are shown in dashed black lines.n = 58 bones from 22 RA participants.(C) Regression plot comparing cortical interruption volume between implementations.(D) Bland-Altman plot comparing in cortical interruption volume between implementations.n = 105 interruptions from 22 participants.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Visualization of cortical break interruption volume and erosion volume in an XCTII image of the 2nd metacarpal.

6 .
Regression plot (A) comparing the erosion volume assessed with MIAF and BAM-Erosion.Bland-Altman plot (B) comparing the differences in erosion volume assessed with MIAF and BAM-Erosion.The mean difference shown as a solid black line and 95 % limits of agreement are shown in dashed black lines.n = 39 erosions from 17 RA participants which came from 5 P2, 25 MC2, 3 P3, and 6 MC4 bones.