Segmentation of Articular Cartilage and Early Osteoarthritis based on the Fuzzy Soft Thresholding Approach Driven by Modiﬁed Evolutionary ABC Optimization and Local Statistical Aggregation

: Articular cartilage assessment, with the aim of the cartilage loss identiﬁcation, is a crucial task for the clinical practice of orthopedics. Conventional software (SW) instruments allow for just a visualization of the knee structure, without post processing, o ﬀ ering objective cartilage modeling. In this paper, we propose the multiregional segmentation method, having ambitions to bring a mathematical model reﬂecting the physiological cartilage morphological structure and spots, corresponding with the early cartilage loss, which is poorly recognizable by the naked eye from magnetic resonance imaging (MRI). The proposed segmentation model is composed from two pixel’s classiﬁcation parts. Firstly, the image histogram is decomposed by using a sequence of the triangular fuzzy membership functions, when their localization is driven by the modiﬁed artiﬁcial bee colony (ABC) optimization algorithm, utilizing a random sequence of considered solutions based on the real cartilage features. In the second part of the segmentation model, the original pixel’s membership in a respective segmentation class may be modiﬁed by using the local statistical aggregation, taking into account the spatial relationships regarding adjacent pixels. By this way, the image noise and artefacts, which are commonly presented in the MR images, may be identiﬁed and eliminated. This fact makes the model robust and sensitive with regards to distorting signals. We analyzed the proposed model on the 2D spatial MR image records. We show di ﬀ erent MR clinical cases for the articular cartilage segmentation, with identiﬁcation of the cartilage loss. In the ﬁnal part of the analysis, we compared our model performance against the selected conventional methods in application on the MR image records being corrupted by additive image noise.


Introduction
The articular cartilage is focused on a resistance to the compressive forces, load distribution, and with a combination of the synovial fluid, frictionless movement of the articular joint components. Structurally, the articular cartilage contains approximately 70-80% of fluid and 20-30% of cellular matrix of the chondrocytes, having a sparse distribution. The chondrocytes are responsible for homeostatic and repair processes, modulating composition of the fluid-like macromolecular network [1,2].
From the view of the clinical practice, osteoarthritis (OA) represents one of the most prevalent musculoskeletal disorders, causing a substantial part of the elderly population. Furthermore, there is an unfavorable prediction of the next increase in the future time due to the population aging [1][2][3]. A severe complication, which relate to OA are structural changes of the articular cartilage, causing its degradation, having progressive and irreversible consequences. Also, it is closely linked with adjacent joint structures pathologies, including, for instance, the subchondral bone and meniscus [4][5][6].
In the clinical diagnosis of the OA, the structural diagnosis is essential. This assessment is based on the indication of a definite osteophyte in anterior-posterior X-ray images. Such image modalities are capable of finding the joint space, but they are not sensitive to a progression of the cartilage loss. This method is frequently used for indication of the effectiveness of disease-modifying OA drugs [7][8][9][10]. Another imaging alternative for the articular cartilage is the ultrasound. [11] As well as the radiographs, ultrasound imaging is not capable of reliably imaging the articular cartilage morphological structure. Thus, tiny structural changes cannot be properly investigated [12,13].
In comparison with others, magnetic resonance imaging (MRI) enables a structural visualization of all the tissues related with the OA disease, including the articular cartilage. A great benefit of MRI is qualitative and quantitative morphologic assessment [14,15]. MRI-based approaches enable characterization and quantification of the articular cartilage biochemical composition. Such techniques include the relaxometry measurements, including T2, T1, T2* imaging and T1rho mapping, sodium imaging, delayed gadolinium-enhanced MRI of cartilage (dGEMRIC), glycosaminoglycan specific chemical exchange saturation transfer (gagCEST), and diffusion weighted imaging (DWI). Compositional MR-based imaging techniques have a strong potential to serve as a quantitative, reproducible, and non-invasive techniques for OA investigation [16][17][18]. On the other hand, current needs of articular cartilage clinical imaging are focused on autonomous extraction and modeling of the cartilage structure, with the regard of the OA. Such segmentation techniques should be capable of extracting the physiological structure of the articular cartilage, and classify the pathological tissues affected by OA. A summary drawback of the MR imaging is a relative smaller proportion of the articular cartilage, regarding the whole knee area. Therefore, the early pathological cartilage loss is insufficiently differentiated from the healthy cartilage, as it is shown in the Figure 1. degradation, having progressive and irreversible consequences. Also, it is closely linked with adjacent joint structures pathologies, including, for instance, the subchondral bone and meniscus [4][5][6].
In the clinical diagnosis of the OA, the structural diagnosis is essential. This assessment is based on the indication of a definite osteophyte in anterior-posterior X-ray images. Such image modalities are capable of finding the joint space, but they are not sensitive to a progression of the cartilage loss.
This method is frequently used for indication of the effectiveness of disease-modifying OA drugs [7][8][9][10]. Another imaging alternative for the articular cartilage is the ultrasound. [11] As well as the radiographs, ultrasound imaging is not capable of reliably imaging the articular cartilage morphological structure. Thus, tiny structural changes cannot be properly investigated [12,13].
In comparison with others, magnetic resonance imaging (MRI) enables a structural visualization of all the tissues related with the OA disease, including the articular cartilage. A great benefit of MRI is qualitative and quantitative morphologic assessment [14,15]. MRI-based approaches enable characterization and quantification of the articular cartilage biochemical composition. Such techniques include the relaxometry measurements, including T2, T1, T2 * imaging and T1rho mapping, sodium imaging, delayed gadolinium-enhanced MRI of cartilage (dGEMRIC), glycosaminoglycan specific chemical exchange saturation transfer (gagCEST), and diffusion weighted imaging (DWI). Compositional MR-based imaging techniques have a strong potential to serve as a quantitative, reproducible, and non-invasive techniques for OA investigation [16][17][18]. On the other hand, current needs of articular cartilage clinical imaging are focused on autonomous extraction and modeling of the cartilage structure, with the regard of the OA. Such segmentation techniques should be capable of extracting the physiological structure of the articular cartilage, and classify the pathological tissues affected by OA. A summary drawback of the MR imaging is a relative smaller proportion of the articular cartilage, regarding the whole knee area. Therefore, the early pathological cartilage loss is insufficiently differentiated from the healthy cartilage, as it is shown in the Figure 1. In this paper, we are focused on the regional segmentation of the knee MR images. Particularly, we are focusing on extraction and modeling of the physiological structure of the articular cartilage, with regard of identification of the early cartilage loss caused by OA. We propose a multiregional segmentation scheme, which is based on soft fuzzy thresholding, where forming the individual segmentation classes is driven by the modified evolutionary artificial bee colony (ABC) optimization scheme. This hybrid segmentation method, including fuzzy based thresholding and evolutionary computing, offers promising results of the cartilage segmentation even in a noisy environment. We provide testing and evaluation of the proposed segmentation method, applied on the native clinical Symmetry 2019, 11,861 3 of 23 MR knee images, containing OA, as well as the synthetically degraded MR images by deterministic noise generators to evaluate robustness of the proposed method.
The organization of the paper is as follows. In Section 2 we describe the recent methods for the cartilage segmentation, which are grouped based on the level of a user's interaction. Section 3 deals with the proposal of the segmentation model for the cartilage segmentation. We define the basic concept of the soft thresholding, including intensity segmentation and spatial aggregation. Further, we introduce the optimization procedure of the soft thresholding based on the modified ABC algorithm. Section 4 brings the results of the proposed model for the cartilage segmentation. We are especially focused on the healthy cartilage segmentation and indication of the early cartilage loss, caused by the early osteoarthritis. Section 5 deals with the quantitative comparison and objective evaluation of the proposed method against selected state of the art segmentation techniques. Section 6 brings the discussion about the results of the proposed segmentation scheme for the cartilage segmentation.

Recent Methods of Cartilage Segmentation
This chapter describes recent segmentation models for the cartilage segmentation. Generally, there are various reasons for the cartilage segmentation, as it is modeling of the morphological cartilage structure, or the pathological finding objective assessment. On the one hand, the ideal segmentation method should work as autonomous as possible, on the other hand, different MR sequences provide cartilage imaging with different intensity manifestation. Therefore, a development of the fully autonomous system allowing for precise segmentation is still the topic of the recent development. Individual methods are grouped by level from the user's interaction [19,20].
The manual segmentation, performed by the clinical experts, has been perceived as a gold standard for the segmentation method's performance assessment. It is a substantially demanding and time-consuming approach which may represent analysis over several hours depending on the image size, quality, and particular pathological findings. Furthermore, this methods lack of relevant reproducibility is caused by subjective error and requires special training [21,22].
From a practical point of view, a balance between the algorithm's performance, preciseness, robustness, and the user interaction level should be kept. An ideal automatic segmentation algorithm would restrict the user's interaction to a minimum and, at the same time, maximize the algorithm's preciseness. There are several factors making this task complicated. The articular cartilage thickness is thin, less than a millimeter. Also, a significant issue is having sufficient contrast achieved between the cartilage and adjacent tissues. Furthermore, the cartilage tissue is not homogenous in certain parts. Therefore, the automatic algorithm should be trained in order to recognize these patterns. The fully cartilage segmentation includes the following techniques: Texture analysis [22], supervised learning [23], statistical methods (active shape models and adaptive template matching) [21], graph-cut methods [22,24], and edge detection methods [19,20,25,26].
The semi-automatic segmentation methods are aimed to reduce the user's interaction as much as possible. These methods commonly require the initialization, meaning that the user is required to specify certain parameters driving the segmentation process. The semi-automatic methods may be divided as follows: Intensity-based [27], thresholding, watershed [28], edge detection [29], energy minimization [30], Live Wire [31], and active contours [30,32].

Materials and Methods
In the comparison with the adjacent knee tissues, the articular cartilage manifestation belongs to a different part of the intensity range. By this assumption, it can be modelled by the regional segmentation identifying individual tissues by their intensity spectrum. Since different MR devices and sequences produce the knee images with various resolutions and image quality, there is not a unified finite intensity range reliably approximating the cartilage area. Furthermore, we must consider that the presence of the image noise and the MR artifacts will significantly modify the intensity distribution, which can never be excluded from the MR image records. The proposed segmentation Symmetry 2019, 11, 861 4 of 23 method compensates these facts by the way that creates adjustable image regions on the base of the pixel's intensity and spatial information driven by evolutionary optimization with the cartilage image features. The overall flow chart of the proposed method is depicted in Figure 2.
to a different part of the intensity range. By this assumption, it can be modelled by the regional segmentation identifying individual tissues by their intensity spectrum. Since different MR devices and sequences produce the knee images with various resolutions and image quality, there is not a unified finite intensity range reliably approximating the cartilage area. Furthermore, we must consider that the presence of the image noise and the MR artifacts will significantly modify the intensity distribution, which can never be excluded from the MR image records. The proposed segmentation method compensates these facts by the way that creates adjustable image regions on the base of the pixel's intensity and spatial information driven by evolutionary optimization with the cartilage image features. The overall flow chart of the proposed method is depicted in Figure 2. In the proposed method, we tried to avoid using the image preprocessing operations to not modify native clinical information. We only use a selection of region of the interest (RoI) in order to In the proposed method, we tried to avoid using the image preprocessing operations to not modify native clinical information. We only use a selection of region of the interest (RoI) in order to expand the area of the articular cartilage due to its small size when comparing with the other knee structures like bones. Unlike many methods for the cartilage segmentation proposed in the recent literature dealing with the cartilage extraction, the proposed method is focused on the cartilage surface analysis, which enables the mapping of tiny structural changes, representing the early cartilage loss.
The basic version of the soft thresholding in [33] utilizes a decomposition of the image area based on the fuzzy sets with clustering and consequent statistical aggregation. This approach represents a general regional image segmentation procedure, not regarding any particular tissues. We bring an optimization of the soft thresholding based on evolutionary computing, selecting the best configuration of the fuzzy class distribution. It is a hybrid segmentation scheme, utilizing the soft thresholding approach which is driven by the evolutionary optimization. Furthermore, in the proposed approach, the configuration of the fuzzy sets are not only done by the clustering, but the real cartilage features are incorporated. This improves accuracy and robustness of the proposed segmentation model.

Soft Segmentation of Intensity Spectrum
A soft multiregional segmentation supposes that the knee image area may be decomposed into a predefined number of the segmentation classes, which should correspond with individual tissues, having different intensity spectrum. In order to perform this task, the MR image histogram is approximated by a sequence of the triangular fuzzy membership functions. In this configuration, the image histogram values are approximated by the fuzzy membership values, so that each pixel will have a certain membership in each segmentation class. We process the image data, represented by 256 gray levels, so the fuzzy approximation is done by: (0; 255) → (0; 1) . Membership values (0; 1) assign a level of the membership for each pixel from the range (0; 255) of gray levels. In this scheme, membership 0 stands for a pixel that does not belong to this class, contrarily membership 1 indicates the full-assignment to this class. In this scheme, the original pixel's distribution expressed by the image histogram is replaced by the system of membership functions (Figure 3).

Soft Segmentation of Intensity Spectrum
A soft multiregional segmentation supposes that the knee image area may be decomposed into a predefined number of the segmentation classes, which should correspond with individual tissues, having different intensity spectrum. In order to perform this task, the MR image histogram is approximated by a sequence of the triangular fuzzy membership functions. In this configuration, the image histogram values are approximated by the fuzzy membership values, so that each pixel will have a certain membership in each segmentation class. We process the image data, represented by 256 gray levels, so the fuzzy approximation is done by: ( ; ) → ( ; ). Membership values ( ; ) assign a level of the membership for each pixel from the range ( ; ) of gray levels. In this scheme, membership 0 stands for a pixel that does not belong to this class, contrarily membership 1 indicates the full-assignment to this class. In this scheme, the original pixel's distribution expressed by the image histogram is replaced by the system of membership functions (Figure 3). For construction of such a segmentation model, the triangular functions are defined. Individual triangular membership functions are defined based on its centroid. In this case, a centroid is represented by the triangular membership vertex. Supposing we have a system of p triangular membership classes: ( ), = , , … , , we need to determine an equivalent system of centroids represented by the triangular feature vector: = , , … , , where V represents a vertex of a respective triangular function. An illustration of a system of the triangular functions approximating the MR image histogram is depicted in Figure 3.
A membership of pixel k in the image I(k), in the n th segmentation class is given as ( ). When using the triangular function, a membership satisfies a condition: Thus, each pixel is transformed into a space of the fuzzy membership, defined as follows: ( ( )) = ( ( )) ( ( )) … ( ( )) (2) For construction of such a segmentation model, the triangular functions are defined. Individual triangular membership functions are defined based on its centroid. In this case, a centroid is represented by the triangular membership vertex. Supposing we have a system of p triangular membership classes: µ n (x), n = 1, 2, . . . , p, we need to determine an equivalent system of centroids represented by the triangular feature vector: V = V 1 , V 2 , . . . , V p , where V represents a vertex of a respective triangular function. An illustration of a system of the triangular functions approximating the MR image histogram is depicted in Figure 3.
A membership of pixel k in the image I(k), in the n th segmentation class is given as µ n (k). When using the triangular function, a membership satisfies a condition: Thus, each pixel is transformed into a space of the fuzzy membership, defined as follows: µ(I(k)) = µ 1 (I(k)) µ 2 (I(k)) . . . µ p (I(k)) When using the triangular fuzzy membership functions, only two adjacent elements of µ(I(k)) are non-zero. A system of the triangular membership functions satisfies the following requirements: Normality: max(µ n (k)) = 1.

Process of Centroids Extraction
The original idea of the centroid extraction deals with the clustering methods. Nevertheless, the clustering methods require a choice of the initial centroid selection. When selecting improper initial centroids, it may result to inaccurate segmentation results, which badly reflect the knee structures. Note that there is not a versatile method for the initial cluster's selection.
In the proposed method, we use a hybrid scheme for the centroid's definition, based on the K-means clustering and modified ABC evolutionary optimization algorithm. From the K-means, we only use the found centroids, consequently serving for the initial forming of the triangular fuzzy functions, not representing the articular cartilage. The physiological cartilage centroid, as well as others, are found by a modified version of the ABC algorithm utilizing the physiological cartilage features.

Modified ABC Algorithm
ABC algorithm represents an evolutionary optimization algorithm inspired by a bee swarm seeking for the food [34][35][36][37]. In this algorithm, there are three groups of artificial bees: Employed bees (EB), onlooker bees (OB), and scout bees (SB). Each kind of the bees search for food sources. A food source represents one solution of the optimization problem. Each such solution particularly represents one configuration of the fuzzy class's distribution for the image segmentation. A proper distribution of the fuzzy classes is a crucial step in the whole segmentation algorithm. Within the optimization procedure, we are searching for the configuration of the fuzzy sets having the maximal compactness. Therefore, we use the entropy for an evaluation of each fuzzy sets configuration. In our approach, we use the Kapur's entropy as an evaluator of each configuration. In this context, we suppose that a higher entropy corresponds with better fuzzy sets distribution. The optimization problem is determined by finding the optimal distribution of fuzzy classes centroids based on the fitness function. The fitness function is based on the Kapur's entropy calculation for each region. The optimization process searches for maximization of the entropy function of the segmentation model. For this task we use the fitness function as an evaluator of individual solutions.
Firstly, we suppose an equal number of EB and OB. We also suppose that one EB has one food source. The parameter X i = X i,1 , X i,2 , . . . , X i,p represents i th solution in the bee swarm, where p stands for a size of the dimension, and in other words these are the number of optimized centroids of the multiregional soft segmentation model.
In this point, we specify X i solution in the bee's swarm. In the conventional ABC algorithm, these solutions are given randomly. In the proposed method, we utilize the real cartilage features in a combination with the random number generator. In order to determine representative cartilage features, we use the estimator of the physiological cartilage intensity average value. In order to perform this task, we employed the segmentation model based on the active contours driven by minimization of the Gaussian energy [32]. This model can adopt the healthy cartilage structure. Since the healthy cartilage is represented by the nearly focused intensity spectrum, not containing intensity fluctuations, it may be well approximated by the Gaussian distribution. Figure 4 shows an example of the physiological cartilage modeling by using the active contour driven by Gaussian energy minimization.  We analyzed the sample of the MR knee images containing 260 image records. 200 of these images were used for investigation of the healthy cartilage features, the rest of them were used for testing of the proposed model with the target of the osteoarthritis identification. All the images were acquired in MR 1.5T with the fat suppression technique, within the period 2014-2017. Images have the unified spatial resolution 1200 × 800 px and dynamical range ( ; ) gray levels. All the images were acquired in the DICOM format. For each of them we obtained the binary model reflecting the physiological structure of the cartilage. Based on these models, we compute an estimation of the intensity spectrum based on two alternative methods for interval estimation. Firstly, we selected the median interval estimation with the reliability 95%. Estimation of median range for the physiological cartilage is given: We analyzed the sample of the MR knee images containing 260 image records. 200 of these images were used for investigation of the healthy cartilage features, the rest of them were used for testing of the proposed model with the target of the osteoarthritis identification. All the images were acquired in MR 1.5T with the fat suppression technique, within the period 2014-2017. Images have the unified spatial resolution 1200 × 800 px and dynamical range (0; 255) gray levels. All the images were acquired in the DICOM format. For each of them we obtained the binary model reflecting the physiological Symmetry 2019, 11, 861 7 of 23 structure of the cartilage. Based on these models, we compute an estimation of the intensity spectrum based on two alternative methods for interval estimation.
Firstly, we selected the median interval estimation with the reliability 95%. Estimation of median range for the physiological cartilage is given: where x p represents 100 p% quantile and n represents the range of sample data (we analyzed 200 MR image records). The second alternative, for the physiological cartilage average value estimation, is the Gastwirth median. The Gastwirth median is calculated based on the selected median by the following expression: Estimation of the Gastwirth median is given: We computed the following parameters for median estimation: The median interval estimation (Equation (3) Based on the selected interval estimations, we computed a median estimation for the physiological cartilage from 200 MR images (Table 1). Results from both of the methods do not exhibit significant statistical differences. For the ABC algorithm, we use an intersection of both intervals reported in Table 1. Note that we worked with the MR images represented by 256 shade levels. Based on the results from Table 1, we determine estimation of the physiological cartilage as follows: I c = 209; 211 . Firstly, we define a scheme for the X i in the Equation (6).
where R i,L ∈ 0; 1 stands for the random numbers of i th solution and L denotes a vector of segmentation classes: L = 1, 2, . . . , p . Vector C L stands for the i th sequence of the centroids generated by the K-means clustering in the form: In each C L , we have to identify and replace centroid representing the cartilage area. In order to do this task, we use the following c . C cartilage replaces a respective centroid from C L for each solution in the swarm based on the minimal distance where I R i c stands for randomly chosen value from I c . When X i is generated, each OB generates a new candidate solution V i in the X i neighborhood. V i is given as following: where X jk represents a random candidate solution (i j), k represents a random index, given: k ∈ 1, 2, . . . p and function φ ik represents a generator of the random numbers in the range: Based on the selection, the X i and V i are compared on the base of their fitness functions. If fit V i > fit X i , V i is stored in the memory. Otherwise, a new V i is generated. The maximal number of such repetitions is controlled by the choice limit L v (we use: L v =10). When the L v is exhausted, the X i is perceived as an exhausted food source.
In the second stage of the ABC, OB test and evaluate the individual solutions from a global view. This process is done based on the probabilistic selection (Equation (8)), performed as many times, as many solutions we have for X i .
It is supposed that more optimal solution has a greater P i . When the X i is selected, the V i is determined from Equation (7) and their fitness functions are compared. As same as in the first stage, a solution having greater fitness function is selected.
The last part of the algorithm are the scouts, seeking for the new food sources, instead of exhausted food sources. When flagging an exhausted food source as X i,e , consequently scout finds a new one, and the whole procedure is repeated within a predefined number of iterations (NI) (we use NI = 200). A new food source is given by the expression: where lb, ub represent lower and upper limit of dimension, respectively. In the proposed segmentation method these limits represent boarders of the normalized intensity spectrum (0;1). In the final step, a solution having the biggest P i is selected as the most optimal solution for the centroids specification.

Features Extraction Based on Fitness Function
Generally, it is supposed that for each segmentation class theoretically unlimited alternatives may be defined, better or worse reflect features of a respective knee structure. Each alternative should be considered based on the histogram features. This situation (for V k segmentation class) is illustrated on the Figure 5. The key aspect of the entire genetic optimization is a fitness function, which globally evaluate a respective segmentation class in order to determine an optimal placement of the segmentation classes to reliably approximate individual knee tissues.
Generally, it is supposed that for each segmentation class theoretically unlimited alternatives may be defined, better or worse reflect features of a respective knee structure. Each alternative should be considered based on the histogram features. This situation (for segmentation class) is illustrated on the Figure 5. The key aspect of the entire genetic optimization is a fitness function, which globally evaluate a respective segmentation class in order to determine an optimal placement of the segmentation classes to reliably approximate individual knee tissues. An important aspect is an optimization criterion defining a fitness function. Theoretically, a finite set of features characterizing the articular cartilage intensity manifestation may be used. We suppose that the articular cartilage represents a homogenous compact structure on the MR imaging. When intensity fluctuations are presented, they indicate pathological findings. Based on these facts, Kapur's entropy is used for the fitness function design. It is supposed that the higher entropy respective fuzzy function has the more optimal solution represents. Regarding manifestation of the articular cartilage on the MR images, the MR signal of the physiological cartilage is represented by concentrated and compact intensity values, without significant intensity fluctuations. This fact is also proved by the interval estimations reported in Table 1. In such configurations, it can be supposed that a region of the cartilage, containing the intensity values with a high probability will be approximated with a higher entropy. In this optimization procedure, we search for the maximization of the Kapur's entropy of the segmentation model. The higher entropy we obtain for the particular X , the more suitable solution is. Firstly, the probability of each gray level ( ) is represented by the relative frequency divided by a total number of gray levels, defined as follows: We suppose that the MR image contains L intensity components, where: = , , , … , − . Each segmentation class represented by the triangular function is described by two thresholds t and centroid V in the form: , , , , where i stands for order of the segmentation class. The Thresholding vector T of the segmentation model, containing p classes is given: = , , , , , , , , , , … , , , , , An important aspect is an optimization criterion defining a fitness function. Theoretically, a finite set of features characterizing the articular cartilage intensity manifestation may be used. We suppose that the articular cartilage represents a homogenous compact structure on the MR imaging. When intensity fluctuations are presented, they indicate pathological findings. Based on these facts, Kapur's entropy is used for the fitness function design. It is supposed that the higher entropy respective fuzzy function has the more optimal solution represents. Regarding manifestation of the articular cartilage on the MR images, the MR signal of the physiological cartilage is represented by concentrated and compact intensity values, without significant intensity fluctuations. This fact is also proved by the interval estimations reported in Table 1. In such configurations, it can be supposed that a region of the cartilage, containing the intensity values with a high probability will be approximated with a higher entropy. In this optimization procedure, we search for the maximization of the Kapur's entropy of the segmentation model. The higher entropy we obtain for the particular X i , the more suitable solution is. Firstly, the probability p k of each gray level h(k) is represented by the relative frequency divided by a total number of gray levels, defined as follows: We suppose that the MR image contains L intensity components, where: k = 0, 1, 2, . . . , L − 1. Each segmentation class represented by the triangular function is described by two thresholds t and centroid V in the form: t i,1 , V i, t i,2 where i stands for order of the segmentation class. The Thresholding vector T of the segmentation model, containing p classes is given: The Kapur's entropy enables measuring compactness and separability of the segmentation classes is defined by the following equations: By this way, the multiregional thresholding is configured as a multi-dimensional optimization problem. Based on the maximization, the fitness function, assessing each solution, is defined with the target of the entropy maximization (Equation (16)). In Figure 6, we present a structure of the ABC evolutionary optimization for the cartilage segmentation. In order to track an evolution of the fitness function on the NI, we report the convergence characteristics for the fitness function (Figure 7). In this regard, an important parameter, influencing the fitness function is a number of the food sources . In Figure 6, we present a structure of the ABC evolutionary optimization for the cartilage segmentation. In order to track an evolution of the fitness function on the NI, we report the convergence characteristics for the fitness function ( Figure 7). In this regard, an important parameter, influencing the fitness function is a number of the food sources X i . The following Table 2 brings the statistical evaluation of the ABC optimization for six segmentation classes. The testing was done for 60 MR knee images containing OA. We present extract of six images. In the testing, we evaluated the best, worst, median, and standard deviation of the fitness function (Equation 16). Another challenging issue of the segmentation procedure is the time complexity of the segmentation. In this regard, we evaluated the time complexity for different image resolution. From this point of the view, size of the image matrix plays an important role (Table 3). We can interpret the intensity segmentation model of the articular cartilage in the form of a sequence of segmentation classes. The outputs a)-f) in Figure 8 represent the segmentation classes one to six of the segmentation procedure. It is apparent that the hard tissues (bones) are mostly classified in the class one. This sequence of the segmentation models indicates the most significant content of the articular cartilage in the segmentation classes five and six. Nevertheless, classes three and four indicate the edge points, having a significant association to the articular cartilage. Based on these facts, uniqueness of the physiological cartilage classification is limited by its surface inhomogeneities. In this context, the MR signal distribution is not homogenous. The following Table 2 brings the statistical evaluation of the ABC optimization for six segmentation classes. The testing was done for 60 MR knee images containing OA. We present extract of six images. In the testing, we evaluated the best, worst, median, and standard deviation of the fitness function (Equation (16)). Another challenging issue of the segmentation procedure is the time complexity of the segmentation. In this regard, we evaluated the time complexity for different image resolution. From this point of the view, size of the image matrix plays an important role (Table 3). We can interpret the intensity segmentation model of the articular cartilage in the form of a sequence of segmentation classes. The outputs a)-f) in Figure 8 represent the segmentation classes one to six of the segmentation procedure. It is apparent that the hard tissues (bones) are mostly classified in the class one. This sequence of the segmentation models indicates the most significant content of the articular cartilage in the segmentation classes five and six. Nevertheless, classes three and four indicate the edge points, having a significant association to the articular cartilage. Based on these facts, uniqueness of the physiological cartilage classification is limited by its surface inhomogeneities. In this context, the MR signal distribution is not homogenous.

Local Statistical Aggregation
Intensity soft segmentation may lack of robustness, especially when the noisy pixels are present. It is supposed that the pixels representing the image noise have significantly different intensity in regards to pixels belonging to the articular cartilage area. This fact often leads to the incorrect classification of the noisy pixels, and the articular cartilage area contains so called blind spots deteriorating the structure of the cartilage model. In the cartilage model, we use a local statistical aggregation scanning the membership values of all the pixels. By considering local statistical features of each segmentation region, the originally assigned membership values may be modified, and pixels may be re-classified, depending on the statistical features of the surrounding pixels.
Supposing we define for k pixel its neighborhood η(k), for such area we can consider local statistical features. In our model we use the median aggregation defined as follow: In the comparison with the average value, the median is invariant in a presence of outliers. Such pixels often represent the image noise, deteriorating the segmentation results. With the median as 50% quantile can reliably approximate the mean value, even if the noise is present. We suppose that the median represents a robust position parameter being invariant against the distant observations. Robustness ensures an effective using for the spatial segmentation which should not be affected by deviating membership values, with the goal of maximization of the segmentation efficiency, especially when the image noise is presented. The local aggregation is implemented by using the cyclic convolution with the kernel η(k), going through all the segmentation regions. We report an example of using the median aggregation ( Figure 9) for pixel's neighborhood of eight pixels.

Local Statistical Aggregation
Intensity soft segmentation may lack of robustness, especially when the noisy pixels are present. It is supposed that the pixels representing the image noise have significantly different intensity in regards to pixels belonging to the articular cartilage area. This fact often leads to the incorrect classification of the noisy pixels, and the articular cartilage area contains so called blind spots deteriorating the structure of the cartilage model.
In the cartilage model, we use a local statistical aggregation scanning the membership values of all the pixels. By considering local statistical features of each segmentation region, the originally assigned membership values may be modified, and pixels may be re-classified, depending on the statistical features of the surrounding pixels.
Supposing we define for k pixel its neighborhood η(k), for such area we can consider local statistical features. In our model we use the median aggregation defined as follow: In the comparison with the average value, the median is invariant in a presence of outliers. Such pixels often represent the image noise, deteriorating the segmentation results. With the median as 50% quantile can reliably approximate the mean value, even if the noise is present. We suppose that the median represents a robust position parameter being invariant against the distant observations. Robustness ensures an effective using for the spatial segmentation which should not be affected by deviating membership values, with the goal of maximization of the segmentation efficiency, especially when the image noise is presented. The local aggregation is implemented by using the cyclic convolution with the kernel η(k), going through all the segmentation regions. We report an example of using the median aggregation ( Figure 9) for pixel's neighborhood of eight pixels. It is obvious when not supposing a presence of the noisy pixels, all the pixels will have great membership values for respective region (Figure 9 a). When noise is presented by the centered pixel, indicated as red (Figure 9 b), such pixel has a significantly lower membership value for a respective region, when comparing with surrounding pixels, due to its different intensity features which can be expected in the image noise. The median local aggregation has the task to modify the originally assigned membership in a dependence of the pixel's surrounding. If the analyzed pixels have similar intensity values, the median aggregation does not rapidly change them, contrarily the noise pixels membership may be substantially modified when considering their surroundings. Such a situation is illustrated in Figure 10 where we report an application of the median aggregation for the situation from Figure 9.

Results
The main intention of the proposed model is a proper localization of the healthy articular cartilage, with regard to the early stage of the cartilage loss caused by the osteoarthritis. MR imaging enables a reliable and contrast visualization of the healthy cartilage nevertheless the early pathological changes are often difficult to recognize by naked eyes.
We did the cartilage modeling on two issues. The first mentioned alternative deals with a selection of the region of the interest where the segmentation procedure is consequently applied. A certain disadvantage of this approach is a spatial limitation of the articular cartilage when only a part of the cartilage is subjected to the segmentation analysis. The second alternative deals with the analysis of the whole knee area where the cartilage is relatively small. In Figure 11, we report the example of the segmentation procedure applied on the MR image, containing the early cartilage loss, which is difficult to observe by naked eyes. We marked the analyzed cartilage area as red, and the early cartilage loss as violet in the segmentation model. The proposed segmentation model enables a differentiation between the physiological cartilage structure, indicated as yellow and the early cartilage loss, which is recognizable as the interruption between the healthy cartilage structures. For testing parameters of the ABC: Number of food sources ( = ) and number of iterations (NI=200). It is obvious when not supposing a presence of the noisy pixels, all the pixels will have great membership values for respective region (Figure 9a). When noise is presented by the centered pixel, indicated as red (Figure 9b), such pixel has a significantly lower membership value for a respective region, when comparing with surrounding pixels, due to its different intensity features which can be expected in the image noise. The median local aggregation has the task to modify the originally assigned membership in a dependence of the pixel's surrounding. If the analyzed pixels have similar intensity values, the median aggregation does not rapidly change them, contrarily the noise pixels membership may be substantially modified when considering their surroundings. Such a situation is illustrated in Figure 10 where we report an application of the median aggregation for the situation from Figure 9. It is obvious when not supposing a presence of the noisy pixels, all the pixels will have great membership values for respective region (Figure 9 a). When noise is presented by the centered pixel, indicated as red (Figure 9 b), such pixel has a significantly lower membership value for a respective region, when comparing with surrounding pixels, due to its different intensity features which can be expected in the image noise. The median local aggregation has the task to modify the originally assigned membership in a dependence of the pixel's surrounding. If the analyzed pixels have similar intensity values, the median aggregation does not rapidly change them, contrarily the noise pixels membership may be substantially modified when considering their surroundings. Such a situation is illustrated in Figure 10 where we report an application of the median aggregation for the situation from Figure 9.

Results
The main intention of the proposed model is a proper localization of the healthy articular cartilage, with regard to the early stage of the cartilage loss caused by the osteoarthritis. MR imaging enables a reliable and contrast visualization of the healthy cartilage nevertheless the early pathological changes are often difficult to recognize by naked eyes.
We did the cartilage modeling on two issues. The first mentioned alternative deals with a selection of the region of the interest where the segmentation procedure is consequently applied. A certain disadvantage of this approach is a spatial limitation of the articular cartilage when only a part of the cartilage is subjected to the segmentation analysis. The second alternative deals with the analysis of the whole knee area where the cartilage is relatively small. In Figure 11, we report the example of the segmentation procedure applied on the MR image, containing the early cartilage loss, which is difficult to observe by naked eyes. We marked the analyzed cartilage area as red, and the early cartilage loss as violet in the segmentation model. The proposed segmentation model enables a differentiation between the physiological cartilage structure, indicated as yellow and the early cartilage loss, which is recognizable as the interruption between the healthy cartilage structures. For testing parameters of the ABC: Number of food sources ( = ) and number of iterations (NI=200).

Results
The main intention of the proposed model is a proper localization of the healthy articular cartilage, with regard to the early stage of the cartilage loss caused by the osteoarthritis. MR imaging enables a reliable and contrast visualization of the healthy cartilage nevertheless the early pathological changes are often difficult to recognize by naked eyes.
We did the cartilage modeling on two issues. The first mentioned alternative deals with a selection of the region of the interest where the segmentation procedure is consequently applied. A certain disadvantage of this approach is a spatial limitation of the articular cartilage when only a part of the cartilage is subjected to the segmentation analysis. The second alternative deals with the analysis of the whole knee area where the cartilage is relatively small. In Figure 11, we report the example of the segmentation procedure applied on the MR image, containing the early cartilage loss, which is difficult to observe by naked eyes. We marked the analyzed cartilage area as red, and the early cartilage loss as violet in the segmentation model. The proposed segmentation model enables a differentiation between the physiological cartilage structure, indicated as yellow and the early cartilage loss, which is recognizable as the interruption between the healthy cartilage structures. One of the segmentation procedure implementations is the analysis of time-inverse knee images, which are evaluated for instance for the reason of presence the focal defects (Figure 12 a)), where the focal defect is manifested on the trochlear socket. The red RoI is segmented with eight classes and aggregated by the median window with size 9 × 9 pixels. In the Figure 12b), we present the resulting cartilage model. One of the segmentation procedure implementations is the analysis of time-inverse knee images, which are evaluated for instance for the reason of presence the focal defects (Figure 12a), where the focal defect is manifested on the trochlear socket. The red RoI is segmented with eight classes and aggregated by the median window with size 9 × 9 pixels. In the Figure 12b, we present the resulting cartilage model.
Next, we tested the segmentation procedure in sagittal 2D SE images. Figure 13a represents the articular cartilage with the cartilage defect (arrows). The segmentation procedure differentiate the cartilage area (blue RoI), where other knee structures are suppressed. Next, we tested the segmentation procedure in sagittal 2D SE images. Figure 13a) represents the articular cartilage with the cartilage defect (arrows). The segmentation procedure differentiate the cartilage area (blue RoI), where other knee structures are suppressed. One type of the articular cartilage morphological assessment is FLASH imaging. This technique is also usable for evaluation of time change of the cartilage thickness and volume. The physiological cartilage characterized by high signal intensity in the FLASH imaging. For this reason, the superficial lesions are better observable in deeper contrast. On the other hand, such lesions may be incorrectly jointed with the cartilage defects, for instance fissures. In the following output ( Figure  14), we present the cartilage model, reflecting the complete cartilage loss of medial area with the degenerative lesions. Next, we tested the segmentation procedure in sagittal 2D SE images. Figure 13a) represents the articular cartilage with the cartilage defect (arrows). The segmentation procedure differentiate the cartilage area (blue RoI), where other knee structures are suppressed. One type of the articular cartilage morphological assessment is FLASH imaging. This technique is also usable for evaluation of time change of the cartilage thickness and volume. The physiological cartilage characterized by high signal intensity in the FLASH imaging. For this reason, the superficial lesions are better observable in deeper contrast. On the other hand, such lesions may be incorrectly jointed with the cartilage defects, for instance fissures. In the following output ( Figure  14), we present the cartilage model, reflecting the complete cartilage loss of medial area with the degenerative lesions. One type of the articular cartilage morphological assessment is FLASH imaging. This technique is also usable for evaluation of time change of the cartilage thickness and volume. The physiological cartilage characterized by high signal intensity in the FLASH imaging. For this reason, the superficial lesions are better observable in deeper contrast. On the other hand, such lesions may be incorrectly jointed with the cartilage defects, for instance fissures. In the following output (Figure 14), we present the cartilage model, reflecting the complete cartilage loss of medial area with the degenerative lesions.
The last application, where we tested the proposed segmentation methodology is the MR imaging with variable force fields. These modalities are not conventionally used, but they are mostly used in an area of the clinical research. The main benefit of stronger magnetic fields is in contrast imaging of the articular cartilage. On the other hand, there is worse magnetic susceptibility in tissues, and images are more inclinable to the flow artefacts. Based on the studies, using the 7.0T MR imaging this modality is more effective for imaging of the morphological structure of the articular cartilage. The benefit of this imaging is achieving a better resolution, with a shorter acquisition time, when comparing with the 3.0T systems.
The following results shows application possibilities of the articular cartilage modeling from 3.0T ( Figure 15) and 7.0T MR images ( Figure 16). Analyzed data shows, in both cases, the physiological cartilage, which is represented by a higher level of the MR signal, which leads to a better contrast in comparison with adjacent tissues. The first important aspect in homogeneousness of the cartilage, which is, in the segmentation models, accompanied mostly by unified color spectrum (yellow spectrum). An important aspect of the segmentation model is also a robustness against the artefacts, caused by the magnetic field. The 7.0T image data are affected by the magnetic susceptibility, this phenomenon is indicated by arrows (Figure 16), nevertheless, this artefact is eliminated in the model. The last application, where we tested the proposed segmentation methodology is the MR imaging with variable force fields. These modalities are not conventionally used, but they are mostly used in an area of the clinical research. The main benefit of stronger magnetic fields is in contrast imaging of the articular cartilage. On the other hand, there is worse magnetic susceptibility in tissues, and images are more inclinable to the flow artefacts. Based on the studies, using the 7.0T MR imaging this modality is more effective for imaging of the morphological structure of the articular cartilage. The benefit of this imaging is achieving a better resolution, with a shorter acquisition time, when comparing with the 3.0T systems.
The following results shows application possibilities of the articular cartilage modeling from 3.0T ( Figure 15) and 7.0T MR images ( Figure 16). Analyzed data shows, in both cases, the physiological cartilage, which is represented by a higher level of the MR signal, which leads to a better contrast in comparison with adjacent tissues. The first important aspect in homogeneousness of the cartilage, which is, in the segmentation models, accompanied mostly by unified color spectrum (yellow spectrum). An important aspect of the segmentation model is also a robustness against the artefacts, caused by the magnetic field. The 7.0T image data are affected by the magnetic susceptibility, this phenomenon is indicated by arrows (Figure 16), nevertheless, this artefact is eliminated in the model.   The last application, where we tested the proposed segmentation methodology is the MR imaging with variable force fields. These modalities are not conventionally used, but they are mostly used in an area of the clinical research. The main benefit of stronger magnetic fields is in contrast imaging of the articular cartilage. On the other hand, there is worse magnetic susceptibility in tissues, and images are more inclinable to the flow artefacts. Based on the studies, using the 7.0T MR imaging this modality is more effective for imaging of the morphological structure of the articular cartilage. The benefit of this imaging is achieving a better resolution, with a shorter acquisition time, when comparing with the 3.0T systems.
The following results shows application possibilities of the articular cartilage modeling from 3.0T ( Figure 15) and 7.0T MR images ( Figure 16). Analyzed data shows, in both cases, the physiological cartilage, which is represented by a higher level of the MR signal, which leads to a better contrast in comparison with adjacent tissues. The first important aspect in homogeneousness of the cartilage, which is, in the segmentation models, accompanied mostly by unified color spectrum (yellow spectrum). An important aspect of the segmentation model is also a robustness against the artefacts, caused by the magnetic field. The 7.0T image data are affected by the magnetic susceptibility, this phenomenon is indicated by arrows (Figure 16), nevertheless, this artefact is eliminated in the model.

Quantitative Comparison and Segmentation Performance
This section is dedicated to the quantitative testing and evaluation of the proposed segmentation model. The testing was done based the gold standard cartilage segmentation, defined by the clinical experts from the orthopedics (Figure 17), and against selected regional segmentation methods. In the second part of the testing, we evaluated optimal settings of the local aggregation for

Quantitative Comparison and Segmentation Performance
This section is dedicated to the quantitative testing and evaluation of the proposed segmentation model. The testing was done based the gold standard cartilage segmentation, defined by the clinical experts from the orthopedics (Figure 17), and against selected regional segmentation methods. In the second part of the testing, we evaluated optimal settings of the local aggregation for the cartilage segmentation. The gold standard images for the comparison were done by manual tracing of the cartilage boarders by clinical experts, having the attestation from orthopedics. and . Segmentation Overlapping: overlapping by is defined as follows: where ( , ) is an overlapping between the regions and defined as follows: We used two complementary descriptors: C( → ) and C( → ). We did testing for 100 MR image records where cartilage was extracted. Averaged results of native images (Nat), as well as additive Gaussian (Gauss) and multiplicative Rayleight (Mult) noise of all the tests are summarized in the Table 4. The best result, for each test, are highlighted. For each MR image, the gold standard was generated ( Figure 17) by using the manual segmentation by clinical expert from the orthopedics.  We tested two alternatives of the proposed soft multiregional segmentation, differing in the local aggregation: Median (MedAg) and average (AvAg) aggregation. The greater values of RI and C we give, the better result of test we obtain, while the lower VI is, the better results we obtain. The reported results in the Table 4 show that the proposed method gives the best, and the most robust results, when comparing with other considered methods, even in the noisy environment represented Firstly, we present a comparative analysis against some conventional segmentation methods. For this task, we consider methods based on the fuzzy approach, as well as the methods taking into account the spatial information. These methods are considered as significantly representative for the comparison. A probability of pixel's belonging to respective class, in a frame of spatial restrictions, is defined as spatial probability.
The following scalar metrics are used for the quantitative comparison: Rand Index (RI): Measures a similarity between two segmentation regions. RI compares the compatibility of an assignment between pairs of elements in two regions. The RI formulation is as follows: where N stands for number of pixels, n 11 denotes a number of pairs belonging to the same area C 1 and C 2 , and n 00 stands for a number of pairs in different segmentation classes. RI gives values 0-1 where 0 indicates completely dissimilar regions and 1 stands for exactly same data. Variation Information (VI): measures distance between two segmentations in a sense of their conditional entropy, which is defined: where H (C i ) represents an entropy associated with the C i and I is a mutual information between C 1 and C 2 .
Segmentation Overlapping: overlapping C 1 by C 2 is defined as follows: where O(R 1 , R 2 ) is an overlapping between the regions R 1 and R 2 defined as follows: We used two complementary descriptors: C(S 2 → S 1 ) and C(S 1 → S 2 ). We did testing for 100 MR image records where cartilage was extracted. Averaged results of native images (Nat), as well as additive Gaussian (Gauss) and multiplicative Rayleight (Mult) noise of all the tests are summarized in the Table 4. The best result, for each test, are highlighted. For each MR image, the gold standard was generated ( Figure 17) by using the manual segmentation by clinical expert from the orthopedics. We tested two alternatives of the proposed soft multiregional segmentation, differing in the local aggregation: Median (MedAg) and average (AvAg) aggregation. The greater values of RI and C we give, the better result of test we obtain, while the lower VI is, the better results we obtain. The reported results in the Table 4 show that the proposed method gives the best, and the most robust results, when comparing with other considered methods, even in the noisy environment represented by the Gaussian and multiplicative noise. When analyzing only the proposed method with different aggregators, we mostly get better results for the median aggregation. This fact is predictable regarding the fact that median represents a robust estimator of position.
In the next step, we compared the median aggregation with a different number of the segmentation classes ( Figure 18). When selecting a higher number of the segmentation classes, the knee segmentation model (blue contour) contains tiny pixels clusters, especially in the femoral bone area, as well as the cartilage is severed. Contrarily, when selecting a lower number of clusters, the early cartilage loss is badly detectable in the model, due to a lower sensitivity.
The next important aspect of the testing is the robustness against the additive image noise. We report results of this testing for MR image RoI 500 × 400 pixels. Analysis is done for the additive image noise salt and pepper (Figure 19), additive Gaussian (Figure 20), and multiplicative Rayleight noise ( Figure 21). regarding the fact that median represents a robust estimator of position.
In the next step, we compared the median aggregation with a different number of the segmentation classes ( Figure 18). When selecting a higher number of the segmentation classes, the knee segmentation model (blue contour) contains tiny pixels clusters, especially in the femoral bone area, as well as the cartilage is severed. Contrarily, when selecting a lower number of clusters, the early cartilage loss is badly detectable in the model, due to a lower sensitivity. The next important aspect of the testing is the robustness against the additive image noise. We report results of this testing for MR image RoI 500 × 400 pixels. Analysis is done for the additive image noise salt and pepper (Figure 19), additive Gaussian (Figure 20), and multiplicative Rayleight noise ( Figure 21).  area, as well as the cartilage is severed. Contrarily, when selecting a lower number of clusters, the early cartilage loss is badly detectable in the model, due to a lower sensitivity. The next important aspect of the testing is the robustness against the additive image noise. We report results of this testing for MR image RoI 500 × 400 pixels. Analysis is done for the additive image noise salt and pepper (Figure 19), additive Gaussian (Figure 20), and multiplicative Rayleight noise ( Figure 21).  area, as well as the cartilage is severed. Contrarily, when selecting a lower number of clusters, the early cartilage loss is badly detectable in the model, due to a lower sensitivity. The next important aspect of the testing is the robustness against the additive image noise. We report results of this testing for MR image RoI 500 × 400 pixels. Analysis is done for the additive image noise salt and pepper (Figure 19), additive Gaussian (Figure 20), and multiplicative Rayleight noise ( Figure 21).  The proposed method was tested for deterministic noise reported in Figure 19, Figure 20, and Figure 21 for variable noise settings from softer deterioration of the intensity distribution, up to the cases where the intensity information completely missing. Based on the comparison, the proposed method gives the most robust results for the multiplicative noise. In this case, the shape parameters of the articular cartilage are mostly perceived. Contrarily, it is apparent that the proposed method is the most sensitive when the Gaussian noise is presented. Despite the cartilage contour smoothness being perceived, the main observable fact is the discontinuity of the cartilage structure where the cartilage loss is not nearly recognizable.  The proposed method was tested for deterministic noise reported in Figures 19, 20, and 21 for variable noise settings from softer deterioration of the intensity distribution, up to the cases where the intensity information completely missing. Based on the comparison, the proposed method gives the most robust results for the multiplicative noise. In this case, the shape parameters of the articular cartilage are mostly perceived. Contrarily, it is apparent that the proposed method is the most sensitive when the Gaussian noise is presented. Despite the cartilage contour smoothness being perceived, the main observable fact is the discontinuity of the cartilage structure where the cartilage loss is not nearly recognizable.
The last tested parameter of the segmentation procedure is the size of the median aggregation window. Based on the pixel's surrounding, pixel's membership may be re-considered. Based on the experimental results, the median aggregation appears as the most effective. Nevertheless, the size of the aggregation window has significant impact on the segmentation effectivity. The greater aggregation window is selected, the more pixels in neighborhood are considered. In order to clarify this situation, we report the comparative analysis for different aggregation settings of the aggregation ( Figure 22).  The last tested parameter of the segmentation procedure is the size of the median aggregation window. Based on the pixel's surrounding, pixel's membership may be re-considered. Based on the experimental results, the median aggregation appears as the most effective. Nevertheless, the size of the aggregation window has significant impact on the segmentation effectivity. The greater aggregation window is selected, the more pixels in neighborhood are considered. In order to clarify this situation, we report the comparative analysis for different aggregation settings of the aggregation (Figure 22). The proposed method was tested for deterministic noise reported in Figures 19, 20, and 21 for variable noise settings from softer deterioration of the intensity distribution, up to the cases where the intensity information completely missing. Based on the comparison, the proposed method gives the most robust results for the multiplicative noise. In this case, the shape parameters of the articular cartilage are mostly perceived. Contrarily, it is apparent that the proposed method is the most sensitive when the Gaussian noise is presented. Despite the cartilage contour smoothness being perceived, the main observable fact is the discontinuity of the cartilage structure where the cartilage loss is not nearly recognizable.
The last tested parameter of the segmentation procedure is the size of the median aggregation window. Based on the pixel's surrounding, pixel's membership may be re-considered. Based on the experimental results, the median aggregation appears as the most effective. Nevertheless, the size of the aggregation window has significant impact on the segmentation effectivity. The greater aggregation window is selected, the more pixels in neighborhood are considered. In order to clarify this situation, we report the comparative analysis for different aggregation settings of the aggregation ( Figure 22).  The aggregation procedure gives a different smoothness of the segmentation classes, depending on the aggregation kernel size. When selecting the 3-square kernel (Figure 22a), the segmentation result gives a relatively higher proportion of the noisy clusters in the femoral bone area (blue contour). Contrarily, when selecting the greatest size 15-square of the median kernel ( Figure 22d) these artifacts are nearly eliminated. Nevertheless, when selecting a higher kernel size, the cartilage is reduced. Such a fact underestimates the effect of the segmentation procedure. This fact is well represented in the Figure 23, indicated as red RoI.

Discussion
The articular cartilage diagnosis is one of the key procedures in the clinical practice of the orthopedics. Therefore, a proper evaluation of the morphological cartilage features, with regards to the early pathological changes is important. The MR imaging represents a current standard for the

Discussion
The articular cartilage diagnosis is one of the key procedures in the clinical practice of the orthopedics. Therefore, a proper evaluation of the morphological cartilage features, with regards to the early pathological changes is important. The MR imaging represents a current standard for the cartilage imaging due to a high spatial resolution and optimization of the imaging by using variable MR sequences. Despite all these benefits, the early cartilage loss, predicting a severe development of the cartilage deterioration, is badly recognizable from the native MR images. For this reason, a development of fully automatic method allowing for the identification and quantification of the physiological cartilage, and early cartilage loss is still challenging issue. The proposed segmentation model utilizes a histogram partitioning based on a sequence of the triangular fuzzy sets. Conventional fuzzy thresholding enables a definition of the fuzzy model based on the centroids calculated by the clustering. Clustering usually requires initialization which influences clusters accuracy. This fact is a limitation for a proper distribution of the intensity segmentation. In the proposed approach, we performed the fuzzy thresholding model driven by the evolutionary optimization, based on the modified ABC algorithm, taking advantage the real cartilage features. Optimization method, in comparison with clustering, generates multiple solutions for the fuzzy set placement, where each of them is evaluated by using the fitness function within the segmentation. Furthermore, the proposed method is able to identify the respective class, representing the articular cartilage. This intensity segmentation is consequently completed by the local spatial aggregation.
In the articular cartilage model, we use the median square kernel window which appears, based on the quantitative testing, as the most robust for both the native MR images and images corrupted by deterministic noise.
In the future, we are going to focus on the cartilage features extraction based on the proposed model. There are two challenging issues. Determining cartilage features, including the cartilage thickness and volumetric parameters enable to evaluate condition of the articular cartilage and they also serve as predictors for objectification of the cartilage deterioration. Since the proposed segmentation model is able to classify the physiological cartilage from the early cartilage loss, it would be worth measuring quantification of this impairment. This model feature would be able to objectivize a cartilage damage level.