Automatic Detection of 3D Joint Erosion in Hand HR-pQCT: A Level Set-Based Approach

—Rheumatoid arthritis (RA) is a chronic inﬂam-matory disease, leading to bone erosion in joints and other complications, which severely affects the life quality of patients. To accurately diagnose and monitor the progression of RA, quantitative imaging and analysis tools are desirable. High-resolution peripheral quantitative computed tomography (HR-pQCT) is such a promising tool for monitoring disease progression in RA. However, automatic analysis tools that can provide useful information using HR-pQCT images still need to be developed. In this paper, we propose a level set-based pipeline for automatic detection of bone erosion. It mainly consists of two steps: joint segmentation and erosion detection. For joint segmentation, we propose an initialization-robust probability-based level set framework, in which Gamma distribution is shown to be feasible for modelling the intensity of joint volumes. For erosion detection, we characterize erosion as concave regions of certain shapes and ﬁt them with a parametric el-lipsoid model. Experiments show that the proposed method outperforms several existing segmentation algorithms and its detection result agrees well with manual annotation.

methods mainly consists of two parts: joint segmentation and erosion detection on the joint contour.The whole procedure is mainly based on multi-view 2D slices, in which erosion detection is mainly performed slice-by-slice in the coronal view with additional information from two other anatomical planes (i.e.transverse and sagittal planes).Only using 2D-projection views for diagnosis is very time-consuming and tedious for radiologists, and important details may be neglected because erosions lie on the 3D joint surface with certain orientations.It is therefore desirable to develop the automatic 3D analyzing tool for erosion detection.
Existing automatic analysis algorithms for joint erosion detection include [9]- [11].Lang et al. [9] follow the standard "segmentation-to-detection" pipeline and make each step fully automatic.For the segmentation part, authors of [9] use a trained model as the initialization to foster the convergence of the 2D active contour model (ACM) [12].For the erosion detection part, with the extracted boundary curve, they model the erosion detection as a classification problem where they use a learned Adaboost classifier to classify points along the curve into erosion versus non-erosion.The other two methods [10], [11] also follow this pipeline with some modifications.Murakami et al. [10] use a deep convolution neural network (CNN) as the classifier.Ren et al. [11] use the level set-based ACM and manual refinement to generate labels for training the segmentation CNN and the Siamese Network to alleviate category imbalance for training the classification models.
Despite the promising performance, the above methods still suffer from the following three issues: • All these methods are based on 2D coronal slice analysis, where intrinsic 3D spatial information is not fully utilized; • In terms of joint segmentation, manual or semi-automatic methods usually choose classical thresholding methods suggested by manufacturer [13], which can give a rough segmentation without manual initialization.Unfortunately, the thresholding methods perform poorly against noises.It has been widely accepted that level set-based ACMs [14]- [17] are good at preserving the topological details during segmentation and robust to noises [18].When they are used in the mentioned automatic detection methods, their refined segmentation results heavily depend on the initialization and usually additional manual annotations are necessary.This is impractical in clinic because it is impossible for radiologists to tediously script 3D initialization case by case; • In terms of erosion detection, the automatic algorithms heavily depend on the labeled training data to learn the well-performing classifier.However, it is challenging to collect enough well-annotated 3D volumes as training samples to customize data-driven models in clinical practice.In this paper, based on the level set methodology [14], [16]- [26], we propose a fully automatic erosion detection algorithm that aims to provide feasible solutions to the issues mentioned in the existing methods.Similar to [9], our method also includes joint segmentation and erosion detection, but the difference is that our method directly processes 3D HR-pQCT data, achieves robust performance against initialization variation in ACM segmentation, and designs effective geometrybased features for erosion detection which does not need extra annotated data to develop data-driven models.Accordingly, our main contribution is threefold: • For joint segmentation, we expound the relationship between the threshold segmentation methods and level set based ones, and a novel initialization-robust level set segmentation framework incorporating probability priors is proposed to achieve stable automatic segmentation of 3D joint volumes; • For erosion detection, we define the erosions on the bone surface as the significant concave regions and voxel-wise surface curvature information is used as features to select candidate erosion voxels.Then we propose to cluster candidate erosion voxels into different regions by their spatial adjacency and use ellipsoid representation as the shape descriptor to conduct region-wise erosion detection; • To our knowledge, our method is the first 3D fully automatic erosion detection algorithms, and its performance can reach consensus with human specialist.The rest of the paper is organized as follows: Section II presents a detailed description of the proposed algorithm.Section III shows the relevant experimental comparison and pathological analysis.Finally Section IV concludes the paper with some discussions.

II. THE PROPOSED MODELS
Fig. 1 shows the workflow of our proposed analysis method.It consists of two parts: 3D joint volume segmentation (II-A) and erosion detection and quantitative analysis (II-B).

A. 3D Joint Segmentation
In this subsection, we will explain the proposed joint segmentation methods in the following order: (and region R 2 ) independently and identically follows the distribution with probability density function (PDF) p(I|R 1 ; θ 1 ) (and p(I|R 2 ; θ 2 )), where θ 1 (and θ 2 ) are unknown distribution parameter vectors.In level set segmentation, the sign of the level set function is used to indicate which part the voxel at x belongs to: Existing level set segmentation methods are categorized into two classes: region-based ones [14], [16], [17] and contourbased ones [19], [20]; and some hybrid works use the two kinds of energy functional together [21], [26].Here we consider the hybrid form, which consists of statistical regional term E R in [17], classical geodesic active contour (GAC) term E GAC in [19], and distance regularization term E reg in [20].The whole energy functional E 0 to be minimized reads where η 1 , η 2 and η 3 are hyper-parameters, G σ is the Gaussian blurring kernel, and the H(•) and δ(•) are the Heaviside function and Dirac function, respectively.H(•) and δ(•) are approximated by the parametric form H (•) and δ (•) mentioned in [19] as follows: When directly applying (2) in 3D joint volume segmentation, the segmentation results are sensitive to initialization.As mentioned before, thresholding methods are still widely used in clinical practice because no specific initialization is needed, but they are sensitive to noises for lack of spatially continuous constraints.Then it is very natural to leverage the thresholding methods to develop an initialization-robust level set formulation.For this purpose, here it is necessary to review the basic idea of global optimal thresholding methods [27].Given the priors of two regions as p(R 1 ) and p(R 2 ), the overall intensity PDF p(I; θ 1 , θ 2 ) reads p(I; θ 1 , θ 2 ) = p(R 1 )p(I|R 1 ; θ 1 ) + p(R 2 )p(I|R 2 ; θ 2 ).(4) Here we assume there is no extra prior which indicates P (R 1 ) = P (R 2 ) = 0.5, and then we can neglect P (R 1 ) and P (R 2 ) terms for simplicity.Normally we cannot analytically determine the global optimal intensity threshold T * , and instead we use the thresholding decision criterion interactively to determine the segmented regions.Concretely, given an initial threshold T (i) for the i-th (i = 0, 1, . . . ) iteration, we first estimate the unknown parameters θ denotes that x should be assigned to R (i+1) 1 when the left side of ( 5) is larger, and to R (i+1) 2 otherwise.Then let us move back to level set framework and assume the current i-th (i = 0, 1, . . . ) level set function is given as 2 , and then the PDF parameter vectors are estimated as θ (i) 1 and θ (i) 2 with maximum likelihood estimation (MLE).If we update the level set function from φ (i) (x) to φ (i+1) (x) from perspective of the thresholding decision criterion (5), then the evolving criterion should satisfy 1) For x with φ (i) (x) ≥ 0 and being assigned to R (i+1) 1 by (5), their level set values φ (i+1) (x) should stay positive; 2) For x with φ (i) (x) ≥ 0 and being assigned to R (i+1) 2 by (5), their level set values φ (i+1) (x) should decrease; 3) For x with φ (i) (x) ≤ 0 and being assigned to R (i+1) 2 by (5), their level set values φ (i+1) (x) should stay negative; 4) For x with φ (i) (x) ≤ 0 and being assigned to R (i+1) 1 by (5), their level set values φ (i+1) (x) should increase.Meanwhile, the level set function should not change infinitely.We can notice from (5) that log p(I(x)|R 2 ) ≥ 0 when φ (i+1) (x) should increase, and vice versa, indicating that level set function φ changes synergistically with the thresholding decision criterion ( 5).Here we use the bounded linear form of ( 5) to update φ as where B(X; M ) = min(|X|, M )sgn(X) bounds each element of X in [−M, M ], and M > 0 is the bounding parameter.For the fixed parameter vectors θ 1 and θ 2 , if the updating term in ( 6) is treated as the gradient descent updating of some energy functional E GR (φ; θ 1 , θ 2 ), i.e., for x ∈ Ω and |φ(x where the second case of ∞ acts as the barrier function to make sure we only consider the solutions with bounded φ. Here we replace the regional energy term E R in ( 2) with E GR to finally get our proposed level set energy functional for initialization-robust joint segmentation as Here we retain the contour-based energy part unchanged because such local information is helpful to explicitly describe boundary.When computing the functional derivative of the regional term E R in (2), the H(•) in (2) will become δ(•), denoting it evolves the isosurface mainly on the contour parts.On the contrary, (7) abandons the Heaviside function term and thus updates the regional term in the whole domain, as shown in (6).In 3D joint segmentation, evolution only near the contour parts is very likely to get stuck if the initialization is not close to the final contour.Formally, Theorem 1 shows that minimizing (7) can produce the same result as E R in (2), while (7) evolves the isosurface among the whole region to achieve robust performance against casual initialization.Theorem 1: Given the fixed distributions p(I(x)|R 1 ; θ 1 ) and p(I(x)|R 2 ; θ 2 ), when φ is bounded by some M < ∞, (7) and E R in (2) share the same minimizer.
Noticeably, in [28], [29], the authors connect the two classical variational image processing techniques, Chan-Vese segmentation models [14] and ROF denoising model [30], and propose the global form of Chan-Vese model, whose regional energy term E gC is defined as where c 1 and c 2 are two constants.In fact, ( 9) is a special case of our formulation if the two likelihood functions of ( 7) are taken as two Gaussian ones with the same variance but different means c 1 and c 2 .Correspondingly, from the framework built for bridging thresholding segmentation and global level set segmentation, ( 9) is analogous to Otsu thresholding methods [27], [31]- [33], while (7) is analogous to more general global optimal thresholding methods (4)- (5).
2) Model joint volume intensity with Gamma distribution: To our knowledge, sophisticated features for precise description of 3D joint volumes are not yet available in the rheumatology community.Here we just use intensity information to model the joint volume.However, intensity-based joint modelling within the probability level set framework faces the following challenges: the joint volumes from HR-pQCT consist of interior marrow, joint bone (sclerotin and periosteum) and exterior muscle tissues.Their boundaries are not clearly visible in HR-pQCT images.In clinical research, the joint bones with marrow inside are taken as a whole.Here we also follow this method and take the bone with inner marrow as the foreground.Taking both capacity and complexity into account on foreground intensity modeling, we desire some representative unimodal distribution.
With 10 well-annotated joint volumes at hand, we inspect the statistical properties of volume intensity with respect to the foreground (joints) and background (tissues).The foreground appears to follow as a positive-skewed distribution.Here we choose several positive-skewed distributions, such as Gamma, Rayleigh, log-normal and Weibull distributions.As none of them can perfectly fit the histogram, we use the Kolmogorov-Smirnov distance [34] to measure the approximation error between a distribution and the foreground histogram.Table I shows that Gamma distribution yields the smallest average Kolmogorov-Smirnov distance for representing the annotated foreground region.Because the histogram of the background region appears symmetric, we just use Gaussian distribution for background modeling.One example is shown in Fig. 2, where we use non-standard Gamma distribution and Gaussian distribution to describe the intensity histogram of foreground and background in our cases.3) Numeral solutions for the segmentation models: Let us instantiate the foreground and background PDF terms in (7) with non-standard Gamma distribution and Gaussian distribution as follows: where m, a and b are the location, shape and scale parameter of the Gamma foreground, respectively; µ and σ 2 are the mean and variance of the Gaussian background.During each iteration, keeping φ fixed, and we estimate (m, a, b) and (µ, σ) of foreground and background with approximate MLE in form of smoothed Heaviside function in (3).With these two groups of parameters, we update the level set function according to (6).As for the Gamma part, we firstly estimate the location parameter m by m := min H (φ(x))>0 I(x) and translate I(x) with m as Ĩ(x) := I(x) − m.MLE of (a, b) of Gamma distribution does not have closed form and needs timeconsuming numerical methods to get an approximate solution.[35] proposes the consistent closed-form approximate MLE estimator of Gamma distribution with the same efficiency as the latent one, and in level set framework it reads where  distribution in level set framework reads The workflow of the initialization-robust segmentation is summarized in Algorithm 1.

4) Post-processing steps:
The raw segmentation volumes need further refinement for surface extraction because of the existence of natural bone vessels and some under-segmented parts, as shown in Fig. 3.A common solution is to apply set theory-based morphological operators to fill the volumes.However, we find that the set theory-based morphology even with moderate template size is easy to cover some small concave regions which could be erosion.To mitigate this problem, we use partial differential equation (PDE) based morphology techniques [36] to preserve more details.Compared with the set theory-based morphology, the PDE-based ones maintain digital scalability and sub-pixel accuracy, and Fig. 4 shows an example demonstrating the detail-preserving advantage of PDE-based morphology in our application.The procedure of the PDE-based close operator is in Algorithm 2.
In clinical practice, radiologists prefer smoother surfaces without changing their original topological properties.To smooth the surface while maintaining its shape, we will minimize the following energy functional where V (x) is the input binary segmentation, φ is initialized as the signed distance map of V , and after several iteration of gradient descent update, the smoothed segmentation is obtained by binarizing the level set result of ( 13).

B. Joint Erosion Detection
In this subsection, our joint erosion detection method will be explained in the following order: • Definition of candidate erosion voxels as locally concave ones based on differential geometry in II-B1; • Isosurface-based curvature estimation within the feasible scale space in II-B2 and II-B3; • Spatial clustering algorithm to aggregate adjacent concave voxels into a whole region in II-B4; • Ellipsoid representation as shape descriptor for the false positive reduction in II-B5. 1) Definition of the candidate erosion voxels based on surface curvature: One can see that erosion areas are concave regions, as shown in Fig. 5.Here we use surface curvature to characterize concave regions [37].Concretely, given a C 2 surface S and an inward positive normal direction, a point is concave if its two principal curvatures k 1 , k 2 , Gaussian curvature K and mean curvature H satisfy K = k 1 k 2 > 0 and H = k1+k2 2 < 0. Furthermore, due to resolution restriction, the curvature radius r 1 = 1 |k1| and r 2 = 1 |k2| should be both larger than 1.In a nutshell, the curvature of the concave candidate erosion points should satisfy 2) Isosurface-based curvature estimation: Given a segmentation volume V and its corresponding surface S, the isosurface φ is derived from the signed distance map of V .To apply differential operators to the digital volume under certain scale priors, we use 3D isotropic Gaussian kernel G σ (x, y, z) = ) to convert the isosurface φ into the C ∞ space as φ σ = φ * G σ , where σ involves the scale priors of the target objects.The curvatures estimation with the gradient and Hessian matrix of φ σ is in the following steps [38]: 1) Given gradient g = (∂ x φ σ , ∂ y φ σ , ∂ z φ σ ) T , compute its normal vector n = − g | g| and P = I − n n T , where I is the identity matrix.
2) Calculate the Hessian matrix H of φ σ as 3) Determination of the feasible scale parameter: The curvature estimation algorithm mentioned above is very sensitive to parameter σ of isotropic Gaussian kernel G σ [38], [40], deciding in which scale space the curvatures are calculated.σ should not be too small to be able to exclude noises and expand isolated points to continuous areas.Meanwhile, σ should be bounded because the most salient erosion parts should be included at least.
Inspired by the scale space representation [41], [42], to determine suitable σ o , φ is smoothed with the monotone increasing scale updated by σ ← √ 2σ, and within each smoothed φ σ , Gaussian and mean curvature are estimated and possible erosion voxels are selected by (14).Because the smoother φ σ is, the fewer parts will be assigned to be bent, and there must be a σ c in which almost no concave point is selected.Because when σ is too small, the proportion of the detected concave points will oscillate and not follow the monotonous change with the scale, and from this perspective, to find σ o we trace back from the σ c and decrease the scale according to σ ← σ √ 2 until the proportion of concave points begin to not increase with the decreasing of scale or the trace-back procedure exceeds the prescribed maximum iteration number.In such a way, the suitable scale parameter can be determined to exclude noises while retaining details.
In our experiments, we implement the above scale searching strategy with the initial scale σ = 0.5 on all test samples, and we find that in most cases samples among the current batch result in σ o = 4 as the feasible scale parameter.Therefore, we just take σ o = 4 as the default scale parameter.When given new batches of samples, it is advisable to carry out such a scale searching procedure to tune the scale parameter.

4) Spatial clustering for distinguishing different erosion areas:
The concave candidate voxels within the scale from II-B3 scatter over several parts, and we need to aggregate them into disjoints regions according to the spatial proximity.After that we can analyze whether each region is erosion or not for the detection purposes.Generally, we do not know the number of concave regions in advance, but the spatial proximity of concave voxels tells that adjacent candidate voxels should belong to the same region.
Here we formulate the clustering algorithm leveraging on the aforementioned spatial proximity as follows: Given a binary erosion indicating a mapping E on domain Ω ⊂ R 3 (x ∈ Ω, E(x) = 1 for erosion), we aim to find a mapping f : In other words, the distance of different erosion regions should be at least r.In clinical diagnose, radiologists usually have empirical knowledge of the distance among erosions.Given a fixed r, we propose a satisfying clustering method in Algorithm 3.

Input: Binary volume E indicating the erosion volume;
minimum distance among different erosions r; Output: Clustering mapping f : Ω → C ⊂ N; cluster number end if 10: end for 11: return f and K. 5) Ellipsoid representation as shape descriptor for false positive reduction: By inspecting the concave regions, we can observe that many natural concave regions on the joint are misclassified as erosion, such as glenoid fossa and two sides of joint bones.Compared with erosion regions, these natural concave regions are much flatter and can be distinguished by their shapes.Moreover, some tiny regions with the similar shape as the erosion ones also should not be classified as true erosion because they are not significant enough in clinical study.Both factors lead to false positives at this raw stage.
Effective false positive reduction is therefore desirable.For the first factor, beyond the voxel-wise curvature feature, priors developed from clinical research for describing the shape of erosion region can be available.Inspired by [4], [5], [43] which use half ellipsoid to represent erosion for identification and quantitative analysis, given a cluster of voxels C indicating candidate erosion, we use the minimal outer ellipsoid of C, also named Löwner-John ellipsoid to represent the erosion as an entirety.Suppose the non-degenerate Löwner-John ellipsoid of 3D point set C parameterized as , where A is some 3 × 3 positive-defined matrix and b ∈ R 3 , we model it as the following convex optimization problem and solve it according to [44]: where det(•) denotes the determinant.Then we use the square root of three eigenvalues of A, √ λ 1 ≥ √ λ 2 ≥ √ λ 3 > 0, the length of three axes of the ellipsoid, as shape indicator to distinguish the true erosion from the false positives, where we treat √ λ 1 as the maximal curvature radium and √ λ 3 as the minimal curvature radium, and their ratio k := λ1 λ3 can serve as a size-invariant measure for shape measurement, acting similarly with the classical concept "shape index" proposed in [45].We adopt the ratio k for false positive reduction, because natural concave structures shared by joints are often extremely flat in ellipsoid representation, corresponding to relatively large k, while the ellipsoid of the true erosion tends to have a relatively small k, indicating those significant locally-bending regions.In our detection system, we set by default that the cluster with k less than 10 should be assigned as erosion.The ablation study on the threshold of k is shown in Table III.For the second factor, the cluster will be excluded if its voxel number is less than the threshold n th .

A. Data Acquisition
The data are acquired in the Department of Medicine & Therapeutics, the Chinese University of Hong Kong, with the approval of the Ethics Committee.The DICOM images are exported from XtremeCT machine (82 µm isotropic resolution).MCP2, MCP3 and MCP4 in each image are cropped into individual volumes for analysis in advance.There are 30 volumes in total and 10 of them are well-annotated by trained radiologists.

B. Experimental Comparisons of Volume Segmentation
The following metrics are chosen to measure the similarity of the delineated ground truth G and segmentation V generated from algorithms: the Dice coefficient D, the false positive (FP) ratio P , the false negative (FN) ratio N and the average Hausdorff distance AH [46], defined as and where d(G, V ) = 1

|G|
g∈G min v∈V ||g − v||.Here the first three metrics are very commonly used.The average Hausdorff distance is adopted because it is more subtle for measuring crisp contour difference and robust to outliers.
When the contour-based and distance regularization terms are fixed, the regional term of the level set segmentation methods for comparison include Chan-Vese model [14], Chan-Vese model in explicit Gaussian distribution [17], Global Chan-Vese model [28], modeling joints with Gaussian and modeling joints with Gamma in the proposed framework, where all the volumes share the same initialization.To be exact, if the size of the volume is h × w × l, we set the initialization as the small cuboids of 0.3h × 0.3w × 0.3l with the same center as the input volume.The quantitative comparison is in Table-II and representative qualitative comparison is in Fig. 6.From Table II, we can see that the FP ratio of the Chan-Vese model and Gaussian Chan-Vese model are relatively large while their FN ratio are very small.This is because these two local level set formulations usually cannot evolve well and their contours do not align with the joint surface, making the error parts mainly consist of false positive ones.Meanwhile, the average Hausdorff distance of these two local ones are much larger, also indicating that the surface of the segmentation volume are not well-aligned, which can be observed in Fig. 6(b)(f).For three level set models in a "global" sense, we can see that these models have a much lower average Hausdorff distance than those local ones in Table II, and from Fig. 6(d (c) Global Chan-Vese model [28] (d) Ours (Gamma foreground) (e) Ground Truth (f) Chan-Vese model [14] (g) Global Chan-Vese model [28] (h) Ours (Gamma foreground) observe that our model with Gamma foreground successfully segments the small concave regions while others fail.

C. Experimental Results of Erosion Detection
In previous related works [9]- [11] of joint erosion detection, performance is measured by pointwise classification precision.In other words, each contour point is only regarded as a classification case.However, there are several drawbacks of sodoing: firstly, because the whole pipeline includes segmentation and erosion detection, this classification precision cannot reveal the errors from the segmentation part; secondly, pointwise classification performance ignores the spatial information and conflicts with the intuition of radiologists.Therefore, it is more reasonable to treat each erosion region as an entirety and check whether it is successfully detected or not.Based on the ellipsoid representation of the candidate erosion regions in II-B5, the average precision (AP) in the object detection task [48] can be adopted to measure the erosion detection performance.For the results V got from ellipsoid axis ratio threshold k, the Interaction Over Union (IOU) ratio with threshold th is adopted to determine whether the erosion G is detected by V : where the successful detection is determined by IOU(V, G) > th.
Here we take IOU threshold th = 0.2 by default.Then the precision and recall are calculated as where TP is the true positive number on the whole test volume set, FP is the false positive number and FN is the false negative number.After that different threshold values are taken and the corresponding Precision-Recall (P-R) values are calculated.Then all the P-R pairs are used to draw the P-R curve, from which mean precision of the model can be estimated.For the P-R pairs with the same recall value, the one with the maximum precision value is left on the figure, then after interpolation, the P-R curve can be got and the area under this curve (AUC) is used to measure the precision of the algorithm.Following the routine above, we sample the P-R values when the ellipsoid axis ratio threshold k = 5, 10, 15, 20, 25, 30, as recorded in Table III.From Table III, we can observe when threshold k = 10 the best balance of precision-recall can be achieved, so we take k = 10 by default.Accordingly, the P-R curve corresponding to Table III is depicted in Fig. 7, where the AUC is 0.793, denoting the probability of the erosion detection results of the algorithm coinciding with the manual annotation [49].Moreover, Fig. 8 shows several erosion detection results in 3D view, where the proposed algorithm successfully finds the general location of the significant erosion throughout the joint surface.

D. Parameter Setting
Besides the scale parameter σ 0 and ellipsoid axis ratio threshold k which have been discussed in the previous II-B3, II-B5 and III-C, here we list the default values of other important parameters.In (8), the time step for updating is set to 1, the coefficients are η 1 = 0.5, η 2 = 6 and η 3 = 0.2 by default, and the maximum iteration number is 16 by default.In Algorithm 2, the numerical step parameter is set to t = 4.In (15), the coefficients are τ 1 = 0.4 and τ 2 = 0.2.In Algorithm 3, the minimum erosion distance is set to r = 4.In II-B5, the second threshold about the voxel number of the region is set to n th = 100.All the volumes are pre-processed to meet the Neumann boundary condition.

IV. CONCLUSION
In this paper, we propose an automatic algorithm for 3D joint erosion in hand HR-pQCT based on the level set The proposed algorithm is fully automatic and interpretable and also provide valuable insights for the description of joint erosions from geometrical perspectives.We plan to experimentally use the current version of the algorithm to aid radiologists with daily clinical trials in the Department of Medicine & Therapeutics, the Chinese University of Hong Kong.At this stage, the size of the detected erosion regions do not exactly coincide with the size of the annotated ones.In the future, we will work on improving the system and make the size measurement more accurate, aiming to develop quantitative analysis methodology for better monitoring the progress of RA.

Proof of Theorem 1
When the parameter θ 1 and θ 2 are fixed, let f (x) := − log p(I(x)|R 1 , θ 1 ) + log p(I(x)|R 2 , θ 2 ).Obviously f (x) is integrable.When φ is bounded by a M > 0, we normalize φ into [0, 1] by ψ ← φ+M 2M and use the normalized level set function ψ to simplify the following deduction, then ( 7) and E R in (2) can be rewritten regardless of the constant multiple M as where 1(•) is the indicator function.Threshold of φ in E R changes from 0 to 1 2 as the threshold of ψ because of the linear normalization.With the coarea formula, E GR can be further rewritten as Taking it into (20), we can get To make the equality in (22) hold, |Σ| should be 0. Therefore, ψ * minimizes E co (ψ; µ) for almost every µ ∈ [0, 1].

•
A new initialization-robust level set-based segmentation framework in II-A1; • Study of the intensity probability distribution of joint volume in II-A2; • Numerical solver for the proposed level set segmentation algorithms in II-A3; • Post-processing for joint surface extraction in II-A4. 1) Initialization-robust probability-based level set segmentation framework: Let I : Ω ⊂ R 3 → R + be a volume intensity function, and volume Ω is segmented by a surface C into two non-overlap regions, joint foreground R 1 and background R 2 .The intensity of each voxel in region R 1

Figure 1 :
Figure 1: Overview of the proposed 3D joint erosion detection system.The upper panel describes the segmentation of the joint volume, consisting of the initialization-robust probability-based level set segmentation framework with Gamma modelling (II-A1, II-A2 and II-A3) and necessary post-processing steps (II-A4).The lower panel describes the erosion detection part, consisting of the primary detection of concave erosion voxels using isosurface-based curvature estimation (II-B1, II-B2 and II-B3), spatial clustering to aggregate concave erosion voxels into erosion region (II-B4), and ellipsoid representation for effective false positive reduction (II-B5).

2 .
Then for each x ∈ Ω, we assign it to the newly-partitioned regions R a, b) are the bias-corrected estimators from the biased ones (â, b), respectively.Parameter estimation of Gaussian

Figure 2 :
Figure 2: Fitting the histogram of foreground with a Gamma distribution (left) and the histogram of foreground with a Gaussian distribution (right).

Figure 3 :
Figure 3: Examples of discontinuity and gaps on the surface of the raw segmentation results: (a) is caused by the natural structure of the bone, while (b) is caused by the errors of the segmentation algorithms.
(a) Discontinuity linking results from set theory-based morphology (b) Discontinuity linking results from PDE-based morphology

Figure 4 :
Figure 4: Comparisons between set theory-based morphology and PDE-based ones for linking the discontinuous parts on the joint surface.PDE-based one is better for preserving important details.

2 F 2 and mean curvature is H = Tr(G) 2 ,
and the shape operator matrix can be written as G = −PHP/| g| 3) Because rank(G) = 2, and the two non-zero eigenvalues are the two principal curvatures, the Gaussian curvature reads K = (Tr(G)) 2 −|G| where Tr is the trace and | • | F is the Frobenius norm.

Figure 5 :
Figure 5: One example shows that the erosion region is clearly characterized by its concavity.Here Paraview [39] is used for visualization.

Figure 6 :
Figure 6: Examples of segmentation results using different level set segmentation methods.The slices are taken from coronal plane.Here ITK-SNAP [47] is used for visualization.

Figure 7 :
Figure 7: P-R curve for measuring erosion detection performance.

Figure 8 :
Figure 8: Two examples of erosion detection results produced from our algorithms, where different erosion regions are in different colors.HereParaview[39] is used for 3D visualization.

Table I :
AVERAGE KOLMOGOROV-SMIRNOV DISTANCE BETWEEN PARAMETRIC DISTRIBUTION AND THE FOREGROUND HISTOGRAM.THE SMALLEST ONE IS UNDERLINED AND BOLD.

Table II :
)(h) we can QUANTITATIVE COMPARISON AMONG REPRESENTATIVE LEVEL SET SEGMENTATION METHODS APPLIED IN 3D JOINT SEGMENTATION, THE SAME POST-PROCESSING PROCEDURE ARE ADOPTED FOR DIFFERENT MODELS.HERE THE METRICS ARE TAKEN AS THE MEAN OF THE TEST SAMPLES, AND FOR DICE COEFFICIENT AND AVERAGE HAUSDORFF DISTANCE THE BEST-PERFORMED ONES ARE UNDERLINED AND BOLD.

Table III :
P-R VALUES FOR THE EROSION DETECTION RESULTS