{'fraction_non_alphanumeric': 0.031435907069948026, 'fraction_numerical': 0.010395140584883212, 'mean_word_length': 4.826340751550529, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiclass histogram-based thresholding using kernel density estimation and scale-space representations

We present a new method for multiclass thresholding of a histogram which is based on the nonparametric Kernel Density (KD) estimation, where the unknown parameters of the KD estimate are defined using the Expectation-Maximization (EM) iterations. The method compares the number of extracted minima of the KD estimate with the number of the requested clusters minus one. If these numbers match, the algorithm returns positions of the minima as the threshold values, otherwise, the method gradually decreases/increases the kernel bandwidth until the numbers match. We verify the method using synthetic histograms with known threshold values and using the histogram of real X-ray computed tomography images. After thresholding of the real histogram, we estimated the porosity of the sample and compare it with the direct experimental measurements. The comparison shows the meaningfulness of the thresholding.


Introduction
A histogram-based thresholding algorithm to perform classification tasks has several advantages compared to a direct clustering methods. It is simple to understand and it is computationally efficient once the histogram has been obtained [11]. Histogram-based methods can be classified into two main categories: parametric and non-parametric. Parametric methods, such as the EM algorithm, are based on estimating the parameters of a given model for the histogram distribution. Recently parametric methods using a mixture of generalized Gaussian distributions were proposed [3,4]. Instead, nonparametric methods aim at optimally thresholding the histogram according to a predetermined criterion [11]. Nonparametric approaches have proven to be robust and more accurate than parametric ones in a number of instances, such as high dimensional data and when clusters deviate substantially from Gaussian distributions [26,20]. Nonparametric approaches include the K-means clustering algorithm and the method developed by Kittler & Illingworth [17]. The authors proposed to threshold the histogram by minimizing the Bayes classification error under the assumption of a bimodal Gaussian distribution of the histogram.
Later, their approach was generalized to a mixture of not-well-separated distributions [7]. A comparative analysis of various nonparametric methods can be found in [24,11].
As the number of thresholds increases, the optimization problem solved in nonparametric methods can become computationally intensive and unstable.
Some of these issues can be potentially overcome through other classes of nonparametric methods, e.g. histogram shape-based thresholding. In this regard, Rosin [27] suggested to define the thresholds by finding "corners" in the histogram plot. The shape-based method requires the threshold values to be clearly reflected in the histogram's shape. Yet, noise and bins width variations may affect the accuracy and stability of shape-based thresholding algorithms. To overcome these issues, various filters can be applied to smooth the histogram.
Tsai [28] suggested to find the maximum curvature of the histogram plot after applying a Gaussian filter. In the method developed by Gilles & Heal [10], the optimal thresholds are found by tracking local minima through a scale-space representation. This approach is based on the observation that threshold values correspond to consistent local minima in the histogram through a scale-space representation. Thresholding along the local minima inherits the advantages of parametric methods but it is computationally more stable and efficient. The local minima thresholding minimizes the Bayes classifier error without explicitly identifying the distribution parameters since it only requires that the clusters are separated by meaningful local minima. However, this may not always be the case, specifically when statistical distributions are not-well-separated and merge into unimodal histograms.
Highly-overlapped classes are frequent in low-contrast image thresholding problems [18], particularly when images produced by gradient-based edge detection algorithms [23], medical X-ray computed tomography (XCT) [16], computer vision [13,6], security and surveillance at low-illumination and low-exposure [21,22,14]. Popular edge detection methods, based on the thresholding of gradient maps, i.e., segmenting the gradient map histogram into two clusters, produce images with low contrast corresponding to unimodal histograms with long tails [23]. The clusterability of highly-overlapped clusters is a subject of discussions [12]. Recent work states that highly-overlapped clusters which fail multimodality test should not be clustered [1].
Nevertheless, histogram-based methods have been proposed to threshold unimodal histograms on two clusters. For example, the thresholding can be achieved by minimizing piecewise linear regression by finding the two segments that fit the descending slope of the histogram [8]. Another method is based on assigning a 2D-point in a receiver operating characteristic (ROC) space to each possible threshold of the histogram without a prior reference image. The optimal point and the required corresponding threshold is then determined in the ROC graph [23]. Although these methods are efficient for binary segmentation of unimodal distributions, they are difficult to extend to the multicluster thresholding case. We refer the reader to the review papers [2,27] on nonparametric thresholding of unimodal distributions for more details. In this work, we propose a new multiclass thresholding method which thresholds along the local minima and is able to handle highly-overlapped distributions.
In Section 2, we provide some generalities regarding histogram modeling and we review scale-space representations and their basic properties. In Section 3, we introduce the three main steps of the thresholding algorithm. In Section 4, we verify the proposed algorithm using a set of synthetic histograms with known threshold values. In Section 5, we threshold low-contrast tomographic images of natural porous sample. The obtained segmentation allows one to estimate the sample porosity, which we compare with experimental measurements. Finally, we summarize our work in Section 6.

Mixture models of histograms
Let us define a histogram H as the set of values given by where h i corresponds to the number of occurrences in the bin interval [t i , t i + ∆t], ∆t being the bin width (where all bins have the same width). A histogram-based segmentation algorithm identifies thresholds to split H into C clusters, where C is a fixed input parameter. To formulate this thresholding problem, we assume that the histogram (1) corresponds to a sampled version of a mixture of K strictly concave unimodal distributions (or modes), H j , j = 1, . . . , K, i.e.
where α j is the a priori probability of occurrence of the j-th component and t is a continuous variable. Several such concave unimodal functions exist in the literature. Here, we will use Gaussian functions, i.e.
For any given pair of consecutive modes {H k , H k+1 }, we define the threshold value τ k by minimizing the Bayes classifier error, i.e τ k is defined such that In the following, we assume the modes H k are well-ordered, i.e their mean values fulfill the condition µ k+1 > µ k where the mean value and the variance of H k are respectively defined by If the modes are well-separated, the threshold value τ k can be easily estimated by finding the position of the local minima (see Figure 1a). Instead, for overlapping modes, the intersection points cannot be resolved by finding local   minima (see Figure 1b) and a different strategy must be adopted.
In Section 3, we present the algorithm which includes the reconstruction of a mixture of Gaussians, the detection of local minima and the identification of the appropriate thresholds.

Scale-space representation
Since our algorithm is based on the notion of scale-space representations, in this section we review how such representations are obtained as well as their basic properties. Given a function f (t), the scale-space representation of f at scale s is given by where * is convolution and G(·, 0, s) are Gaussian functions. As the scale s increases, this process filters lower frequency harmonics of f . As a result, no local extrema can appear for a large enough s. This property is at the core of the method to detect consistent local minima by tracking the "long-life" presence of minima through the scales, as described in [10].

Thresholding algorithm
The algorithm is composed of two steps. In Section 3.

Kernel Density Estimation
First, we approximate the histogram (1) using a kernel density estimation technique corresponding to the minimization problem [5] min where KD is a mixture of Gaussians defined by where the weights B = {β j } N j=1 and the variance σ 2 are unknown parameters to be determined. The variance σ 2 is referred to as the bandwidth of the Gaussian kernel, or as the scale in the scale-space nomenclature.
It is worth emphasizing that we propose to associate one Gaussian to each bin, instead of assuming a small number (generally not known in advance) of

Algorithm 1 Expectation-Maximization Kernel Density Estimation algorithm
Initialize: Gaussians. This choice allows one to the estimation of the number of Gaussians in the mixture. Furthermore, since we assume that all bins have the same width, we use the same variance for each Gaussian in the mixture. This has the advantage of avoiding the collapsing problem in the Expectation-Maximization (EM) algorithm described hereafter.
The minimization problem (8) is known to be ill-posed [15] and it usually requires a regularization term. We propose to solve this minimization problem using an EM-type algorithm since it naturally provides a Tikhonov regularization by maximizing the a posteriori probability [5,9,15]. An additional advantage of such an approach is that noise in the histogram can be accounted for. Denoting the Fast Fourier Transform (FFT) and its inverse with F and F −1 , respectively, the EM-KD algorithm is described by Algorithm 1, where all divisions and multiplications must be understood pointwise. It can be shown that Algorithm 1 is equivalent to the Richardson-Lucy deconvolution algorithm [25].
In Fig. 2, we show several iterations from the proposed KD estimation algorithm.

Local minima detection
Given a mixture KD(t, B, σ 2 ), we now need to extract the position of all local minima. These minima can be detected in a semi-analytical fashion. Since an analytic form for H(t) is available, its derivatives can be calculated and their zeros numerically determined. This procedure corresponds to numerically solve the following set of constraints for t: where the analytic form of the derivatives are given by and In the remainder of the paper, the set of local minima will be denoted M and its cardinality by #M.

Thresholding detection algorithm in scale-space representation
In order to detect the expected thresholds, we use a scale-space representation of our histogram. We start by a scale-space representation of our histogram as outlined in Section 3.3.1, and then we present the overall threshold detection algorithm in Section 3.3.2.
The advantage of having a continuous model is twofold: it allows one to easily obtain, on one hand, the scale-space representation of the histogram by choosing ∆σ 2 > 0, and on the other, an "inverse" scale-space representation by choosing ∆σ 2 < 0 (while satisfying the constraint σ 2 + ∆σ 2 > 0). The estimation of an inverse is a key point in dealing with overlapping distributions in the threshold detection algorithm presented in the next section. The magnitude of ∆σ 2 must be chosen less than ∆t, so to ensure enough resolution while minimizing computational costs.

Thresholding detection algorithm
Once the scale-space representation S KD (t, ∆σ 2 ) is constructed, the expected thresholds to segment the histogram into C clusters need to be identified. This procedure is summarized in Algorithm 2.
Find the set M of local minima of KD(t, B, σ 2 ).

Validation on synthetic histograms
We validate the method against a set of 2185 synthetic histograms. The histograms are generated using a mixture of three Cauchy distributions defined by where c j is the a priori probability of the j-th component, and a j and b j are parameters of the Cauchy distribution. We generate a set of 2185 different distributions by uniformly selecting the mixture parameters with −4 < a j < 4, 0.5 < b j < 2 and 0.2 < c j < 0.5.  that the better the initial resolution, the more efficient the algorithm is.

Application to porosity estimation from XCT images
In this section, we present an application of our algorithm to estimate the porosity of some sample through X-ray Computed Tomography (XCT) images. We apply the proposed thresholding algorithm to a stack of a 1000 XCT gray-scale images of a natural porous sample (the algorithm is applied on a unique histogram obtained from all pixels of all images in the stack). An example of such an image is given in Fig. 5a. All details about the sample as well as the XCT imaging technique are provided in [19].  The pixel intensity of an XCT image is linearly related to the pixel porosity through the map [19] where t is the pixel intensity, t s and t v are two reference points that have to be determined. The reference point t s corresponds to the pixel intensity of the solid phase (φ = 0), and t v to the pixel intensity of the void phase (φ = 1).
Such reference points are identified by thresholding the histogram. Given the knowledge of the properties of the XCT imaging and the sample, it is known that the combined image intensity histogram consists of three clusters: the first cluster represents the void phase (white regions in Fig. 5a), the second cluster the unresolved porous phase (gray regions in Fig. 5a) and the third the solid phase (dark grey regions in Fig. 5a) [19].
The combined histogram, shown in Fig. 6a, has three overlapping modes and the minima cannot be resolved. We apply the proposed method (with ∆σ 2 = 0.01) to the histogram and obtain the thresholds τ 1 and τ 2 as indicated in Fig. 6b. Given τ 1 and τ 2 , the reference points t v and t s can be determined by taking the average within the first and third cluster, respectively, i.e Next, the gray intensity histogram is transformed into a porosity histogram via (15). The average porosity is therefore calculated by computing the first moment of the porosity histogram. The estimated and experimentally measured porosities are 0.185 and 0.178, respectively, with a relative error of 3.9%. For comparison, we follow the same procedure and use the classic K-means clustering algorithm. The obtained porosity value in this case is 0.27 which corresponds to a relative error of 51.6%.    In this paper, we introduce a new approach to segment histograms with overlapping modes and not sufficiently resolved. Our algorithm exploits the properties of scale-space representations to automatically find the appropriate thresholds, and has only two parameters, the number of expected classes and ∆σ 2 . We first test the proposed method on a large set of simulated synthetic histograms.
Then, we also the algorithm to estimate the porosity of a real porous rock sample from XCT images and show its perfomance compared to the K-means algorithm.
The proposed algorithm has an intuitive formulation of the 1D histogram, but the extension to the multivariate histograms is not straight forward. The The future work includes the extension of the algorithm for a multivariate histogram. This work requires the formulation of a new set of criteria that define the thresholding surfaces. Current investigations are focused on approaches to automatically identify the "optimal" number of classes, and reduce the number of the algorithm parameters to one.