Brain MR image segmentation for tumor detection based on Riesz probability distributions

ABSTRACT

This research introduces a new approach using the Riesz mixture model for medical image segmentation, specifically for diagnosing and treating brain tumors. We developed a novel technique for pixel classification based on the Riesz distribution, which is generated using an extended Bartlett decomposition. Our work is pioneering, as there are no existing studies addressing this issue in the literature. We aim to demonstrate the effectiveness of this distribution for brain image segmentation. We used the Expectation-Maximization algorithm to estimate the mixture parameters. To validate our segmentation algorithm, we conducted a comparative study with a recent method based on the Wishart distribution using Matlab software. Experiments with the Riesz mixture model showed that our method produces more intuitive results with a recognition rate of 94.52%. These results confirm the reliability of our method in detecting tumors using both synthetic and real brain images.

KEYWORDS:

• EM algorithm
• magnetic resonance imaging
• segmentation
• Riesz probability distribution
• Wishart probability distribution

1. Introduction

The medical image is a crucial part of the medical healthcare. It is the technique of creating visual representation of the interior of a body for clinical analysis and medical intervention. From this perspective, different medical imaging technologies are invested as radiography, computed tomography, ultrasound and magnetic resonance (MR) imaging.

To identify the affected area whether tumours, cyst, oedema and other abnormalities, we opted for medical image segmentation methods. Segmentation of image is an intrinsic research field of image processing and pattern recognition. It is a supervised or unsupervised method of partitioning an image into regions. It is a complex task that is affected by numerous aspects like noise. As a matter of fact, various image segmentation approaches have been proposed in literature as the Thresholding method which aims to discriminate pixels or voxels through using their grey levels. These thresholds can be either global or local. They are determined by using the histogram of an image. The segmentation of the image may be either carried out by the application of many individual thresholds or by utilising a multi-Thresholding technique. There are also region based techniques which are designed to locate homogeneous areas in the images. They create disjoint regions by combining neighbourhood pixels with homogeneous characteristics. The approaches contours are characterised with the presence of borders between regions. We have also pixel classification technique which classifies similar objects into a single cluster. It is created by using multiple supervised or unsupervised clusterings, like: Fuzzy C-means clustering and K-means algorithm (Zhou and Yang 2020). In addition, we mention the probabilistic method in which hidden Markov model is introduced to model the tissue. Moreover, mixture models are used to represent the data density of brain image (Bhagwat et al. 2021; Zitouni and Tounsi 2021; Tounsi and Zitouni 2021; Jiao et al. 2022). We can also use the finite mixture model for complex image segmentation problem which corresponds to a convex combination of two or more probability density functions (Zitouni et al. 2018, 2023). It is used in different fields such as finance, handwriting recognition, medicine, and insurance. Among finite mixture models, Gaussian mixture model (GMM) is one of the most used tool for image segmentation (Saygılı 2022; Huang and Gou 2023). It provides better classification in terms of Bayes risk than several popular methods. The paper highlights an unsupervised method (Panić et al. 2022; Pu et al. 2023) of image segmentation based on the Riesz mixture model and the Expectation-Maximization (EM) algorithm for the parameters estimation of the proposed model. Recently, The Riesz mixture model was developed by Kessentini et al. (2020). It is a combination of two or more Riesz probability densities which is generated by extending the famous Bartlett decomposition of Wishart random matrices (Bartlett 1934) to the Bartlett decomposition of Riesz random matrices.

From 2001 until now Riesz distributions and its derivative probabilities have occupied a lofty position in terms of multivariate analysis. Numerous researchers were highly interested in this new family of probability distributions and they provided multiple theoretical aspects which allow to master, characterise and therefore to understand the usefulness of this new family subsequently in real applications. The Riesz probability distribution was introduced in 2001 by Hassairi and Lajmi (Hassairi and Lajmi 2001) on the modern framework of symmetric cones and simple Euclidean Jordan algebra. In 2008, Hassairi et al. (2008) developed a characterisation of the Riesz probability distribution in the cone of positive definite symmetric matrices to make them accessible to readers. The Riesz probability distribution represents a significant generalisation of the Wishart distribution and therefore, of the usual gamma distribution in the real line. Multiple interesting properties of the latter are identified and discussed (Hassairi et al. 2007, 2017; Tounsi and Zine 2017). Within this framework, numerous other derived probability distributions are addressed such as the inverse Riesz probability distributions (Tounsi and Zine 2012), the Dirichlet-Riesz (Tounsi and Zine 2017) and the Extended matrix-variate Beta distributions (Tounsi 2020). However, the Riesz and the inverse Riesz probability distributions displays a prominent property namely the positivity set that is imposed on the covariance matrices of a multivariate normal distribution. For this reason, we can adopt these distributions for the estimation of covariance matrix of a multivariate normal distribution. It is in fact more advantageous to use these distributions than the Wishart and the inverse Wishart distributions because they depends on multivariate shape parameters thus offering a higher degree of freedom. Therefore, the Riesz and inverse Riesz probability distributions can be used in a wide variety of problems such as statistical analysis, economics and physics. Thus, the search for a statistical models and algorithms facilitating the investigation of the Riesz and inverse Riesz in real world problems remains a topic of high interest to statisticians. We recall that in this context, recently Kessentini et al. (2020) have suggested an algorithm generating the Riesz and the inverse Riesz probability distributions based on an extension of the Bartlett decomposition (Bartlett 1934). Within the theoretical framework, many works about the Riesz probability properties are developed. Yet, in real-life applications no work has been recorded in terms of demonstrating the applicability of this distribution (Kessentini et al. 2020). Hence, our research is considered as the first scientific research that stands for the pioneering work aiming at proving the efficiency of the Riesz distribution in real life problems such as segmentation of medical imagery. The latter represents a challenging problem due to the complexity of the images, as well as the absence of models of the anatomy that fully capture the possible deformations in each structure. In addition, to estimate the parameters of the Riesz mixture model, Kessentini et al. (2020) adopted the EM algorithm (Dempster et al. 1977). For real-life applications like segmentation of complex structure, parameters estimation is a basic step when implementing completely data. From this perspective, we suggest in this paper an iterative algorithm based on the EM algorithm. The proposed segmentation algorithm has been validated on synthetic data generated by using the Riesz distribution and has been tested on real medical images. The obtained results proved the high performance and the goodness of fit of the proposed method with Riesz distribution compared to a recent method of segmentation based on the Wishart distribution (Hidot and Saint-Jean 2010). Besides, the finite sample performance was evaluated by calculating such evaluation criteria as: Misclassification Rate $(\mathit{MCR})$ (Zitouni et al. 2023), Peak Signal to Noise Ratio $(\mathit{PSNR})$ (Saladi and Amutha Prabha 2017), Adjusted Rand Index $(\mathit{ARI})$ (Sundqvist et al. 2023), Rand Index $(\mathit{RI})$ (Fouad et al. 2017), Mirkin Metric $(\mathit{MI})$ (Correa-Morris et al. 2022) and Hubert Index $(\mathit{HI})$ (Somashekara and Manjunatha 2014). The different evaluation criteria values suggest the reliability of the proposed distribution. The remaining sections of this paper are organised as follows. In Section 2, some notations and definitions are provided within the general framework of positive definite symmetric matrices. Furthermore, we recalled some basic concepts of the Riesz mixture model and an estimation of the parameters using maximum likelihood method. In Section 2 image segmentation using Riesz mixture model is explored and discussed. Additionally a comparative study with respect to the Wishart distribution is carried out. The proposed algorithm was applied with synthetic and real medical images. The different experiments results demonstrated the high performance of the proposed unsupervised algorithm. Moreover, the importance of using Riesz models in real life applications is highlighted. The last section wraps up the conclusion and provides new perspectives for future works.

2. Notations and definitions

To clarify the results of this paper, we first recall some notations and review some characteristic properties concerning the Riesz distributions on the modern framework of positive definite symmetric matrices.

2.1. The Riesz probability distribution

Let $V_r$ be the space of real symmetric $\mathit{r}\times \mathit{r}$ matrices equipped with the scalar product $⟨\mathit{X},\mathit{Y}⟩=\mathrm{tr}(\mathit{XY})$, where $\mathrm{tr}(\mathit{XY})$ is the trace of matrices’ ordinary product. For $\mathit{X}={\left({\mathit{x}}_{\mathit{ij}}\right)}_{1\le \mathit{i},\mathit{j}\le \mathit{r}}$ in $V_r$ and $1\le \mathit{l}\le \mathit{r}$, we provide the sub-matrices

$P_l(\mathit{X})={\left({\mathit{x}}_{\mathit{ij}}\right)}_{1\le \mathit{i},\mathit{j}\le \mathit{l}},$

and

$P_l^{\ast }(\mathit{X})={\left({\mathit{x}}_{\mathit{ij}}\right)}_{\mathit{r}-\mathit{l}+1\le \mathit{i},\mathit{j}\le \mathit{r}}.$

For convenience, $P_l(\mathit{X})$ and $P_l^{\ast }(\mathit{X})$ may be considered a matrices in $V_r$ with zero-filled additional rows and columns.

We denote by ${\mathrm{\Omega }}_{\mathrm{r}}$ the cone of positive definite elements of $V_r$ and by ${\mathcal{T}}_r$ the set of $\mathit{r}\times \mathit{r}$ lower triangular matrices with strictly positive diagonal. Any element $\mathit{Y}$ of ${\mathrm{\Omega }}_{\mathrm{r}}$ admits the Cholesky decomposition $\mathit{Y}=\mathit{T}{T}^{\ast }$, where $\mathit{T}$ is in ${\mathcal{T}}_r$ and $T^{\ast }$ is its transpose. This decomposition is unique, and therefore we can define on $V_r$ the automorphism

$\mathit{\pi }(\mathit{Y})(\mathit{X})=\mathit{TX}{T}^{\ast },$

which may be regarded as the product of $\mathit{X}$ by $\mathit{Y}$, and we define the ‘quotient’ of $\mathit{X}$ by $\mathit{Y}$ as

${\pi }^{-1}(\mathit{Y})(\mathit{X})=T^{-1}\mathit{X}(T^{-1}{)}^{\ast }.$

Definition 2.1. An $\mathit{r}\times \mathit{r}$ positive definite symmetric matrix $\mathrm{X}$ is said to have a matrix-variate Wishart distribution with shape parameter $\mathit{p}>\frac{\mathit{r}-1}{2}$ and scale parameter $\mathrm{\Sigma }$ in ${\mathrm{\Omega }}_{\mathrm{r}}$, if its probability density function is defined by

(1) ${\mathcal{W}}_r(\mathit{p},\mathrm{\Sigma })=\frac{{(det\mathrm{\Sigma })}^{-\mathit{p}}}{2^{\mathit{rp}}{\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{p})}{e}^{-\frac{1}{2}⟨{\mathrm{\Sigma }}^{-1},\mathit{X}⟩}(det\mathit{X}{)}^{\mathit{p}-\frac{\mathit{r}+1}{2}}{1}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{X}),$(1)

where ${\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(.)$ denotes the multivariate gamma function given by

(2) ${\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{p})={\pi }^{\frac{\mathit{r}(\mathit{r}-1)}{4}}\prod _{\mathit{j}=1}^r\mathrm{\Gamma }(\mathrm{p}-\frac{\mathrm{j}-1}{2}).$(2)

The Riesz probability distribution was proposed by Hassairi and Lajmi (2001) as an extension of the Wishart one. It is an interesting probability distribution with outstanding probabilistic properties like independence between blocks (Hassairi et al. 2007), invariance (Hassairi and Lajmi 2001), moments and constancy regression on the mean (Hassairi et al. 2017). Departing from the fact that the Riesz probability distribution verifies the property of set positivity that is imposed on the covariance matrices, the latter can be applied for the estimation of a covariance matrix of a multivariate normal distribution. The definition of the Riesz probability distribution is based on the notion of generalised power in the space of real symmetric $\mathit{r}\times \mathit{r}$ matrices, which reduces to the ordinary determinant in a particular situation.

Definition 2.2. For $\mathit{s}=(s_1,\dots ,s_r)\in {\mathbb{R}}^r$, the generalised power of an element $\mathit{X}\in {\mathrm{\Omega }}_{\mathrm{r}}$ is expressed by

(3) ${\mathrm{\Delta }}_s(\mathit{X})={\mathrm{\Delta }}_1(\mathit{X}{)}^{s_1-s_2}{\mathrm{\Delta }}_2(\mathit{X}{)}^{s_2-s_3}\cdots {\mathrm{\Delta }}_r(\mathit{X}{)}^{s_r},$(3)

where ${\mathrm{\Delta }}_l(\mathit{X})$ is the determinant of the $\mathit{l}\times \mathit{l}$ sub-matrix $P_l(\mathit{X})$ of $\mathit{X}$.

Note that ${\mathrm{\Delta }}_s(\mathit{X})=(det\mathit{X}{)}^p$ if $\mathit{s}=(\mathit{p},\dots ,\mathit{p})$ with $\mathit{p}\in \mathbb{R}$. It is also easy to check that ${\mathrm{\Delta }}_{\mathit{s}+\mathit{s}{\mathit{ }}^{\mathrm{\prime }}}(\mathit{X})={\mathrm{\Delta }}_s(\mathit{X}){\mathrm{\Delta }}_{\mathit{s}{\mathit{ }}^{\mathrm{\prime }}}(\mathit{X}),$ where $\mathit{s}{ }^{\mathrm{\prime }}=(s_1^{\prime },\dots ,s_r^{\prime })\in {\mathbb{R}}^r$. In particular, if $\mathit{m}\in \mathbb{R}$ and $\mathit{s}+\mathit{m}=(s_1+\mathit{m},\dots ,s_r+\mathit{m})$, we have ${\mathrm{\Delta }}_{\mathit{s}+\mathit{m}}(\mathit{X})={\mathrm{\Delta }}_s(\mathit{X})\text{break}(det\mathit{X}{)}^m$. Furthermore, Hassairi et al. (2007) demonstrated that, for all $\mathit{X}$ and $\mathit{Y}$ in ${\mathrm{\Omega }}_{\mathrm{r}}$

(4) ${\mathrm{\Delta }}_s(\mathit{\pi }(\mathit{Y})(\mathit{X}))={\mathrm{\Delta }}_s(\mathit{Y}){\mathrm{\Delta }}_s(\mathit{X}),$(4)

and

(5) ${\mathrm{\Delta }}_s({\pi }^{-1}(\mathit{Y})(\mathit{X}))={\mathrm{\Delta }}_{-\mathit{s}}(\mathit{Y}){\mathrm{\Delta }}_s(\mathit{X}).$(5)

Definition 2.3. An $\mathit{r}\times \mathit{r}$ positive definite symmetric matrix $\mathit{X}$ is said to have a Riesz probability distribution on ${\mathrm{\Omega }}_{\mathrm{r}}$ with shape parameter $\mathit{s}=(s_1,\dots ,s_r)\in {\mathbb{R}}^r$, such that $s_j>\frac{\mathit{j}-1}{2}$ for all $1\le \mathit{j}\le \mathit{r}$ and scale parameter $\mathrm{\Sigma }$ in ${\mathrm{\Omega }}_r$, if its probability density function is indicated by

(6) ${\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })=\frac{1}{2^{\sum _{\mathit{j}=1}^r{s}_j}{\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{s}){\mathrm{\Delta }}_s(\mathrm{\Sigma })}{e}^{-\frac{1}{2}⟨{\mathrm{\Sigma }}^{-1},\mathit{X}⟩}{\mathrm{\Delta }}_{\mathit{s}-\frac{\mathit{r}+1}{2}}(\mathit{X}){\mathbf{1}}_{{\mathrm{\Omega }}_r}(\mathit{X}),$(6)

where

(7) ${\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{s})={\pi }^{\frac{\mathit{r}(\mathit{r}-1)}{4}}\prod _{\mathit{j}=1}^r\mathrm{\Gamma }({\mathrm{s}}_{\mathrm{j}}-\frac{\mathrm{j}-1}{2}).$(7)

The probability distribution ${\mathit{R}}_r(\mathit{s},\mathrm{\Sigma })$ reduces to the Wishart ${\mathcal{W}}_r(\mathit{p},\mathrm{\Sigma })$, when $s_1=s_2=\dots =s_r=\mathit{p}$.

In the following theorem, we recall a significant result concerning the expected value of a Riesz random matrix. For further details, please consult this reference (Hassairi et al. 2017).

Theorem 2.4.

Let $\mathit{X}$ be a Riesz probability with parameters $\mathit{s}=(s_1,\dots ,s_r)\in {\mathit{R}}^r$ such that $s_1\ne 0$ and $\mathrm{\Sigma }\in {\mathrm{\Omega }}_{\mathrm{r}}$. Then, the expected value of $\mathit{X}$ is denoted by

(8) $\mathit{E}(\mathit{X})=-\sum _{\mathit{l}=1}^r(s_{\mathit{r}-\mathit{l}}-s_{\mathit{r}-\mathit{l}+1}){\left(P_l^{\ast }(\mathrm{\Sigma })\right)}^{-1}.$(8)

If $\mathrm{X}$ is a Riesz random matrix, then the distribution of $\mathrm{Y}={\mathrm{X}}^{-1}$ is called the inverse Riesz probability distribution (Tounsi and Zine 2012), which can be regarded as a generalisation of the inverse Wishart.

Definition 2.5.

An $\mathit{r}\times \mathit{r}$ positive definite symmetric matrix $\mathit{Y}$ is said to have an inverse Riesz probability distribution on ${\mathrm{\Omega }}_{\mathrm{r}}$ with shape parameter $\mathit{s}=(s_1,\dots ,s_r)\in {\mathbb{R}}^r$, such that $s_j>\frac{\mathit{j}-1}{2}$ for all $1\le \mathit{j}\le \mathit{r}$ and scale parameter $\mathrm{\Sigma }$ in ${\mathrm{\Omega }}_{\mathrm{r}}$, if its probability density function is determined by

(9) $\mathcal{I}{\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })=\frac{1}{2^{\sum _{\mathit{j}=1}^r{s}_j}{\mathrm{\Gamma }}_{{\mathrm{\Omega }}_r}(\mathit{s}){\mathrm{\Delta }}_s(\mathrm{\Sigma })}{e}^{-\frac{1}{2}⟨{\mathrm{\Sigma }}^{-1},Y^{-1}⟩}{\mathrm{\Delta }}_{\mathit{s}+\frac{\mathit{r}+1}{2}}(Y^{-1}){\mathbf{1}}_{{\mathrm{\Omega }}_r}(\mathit{Y}).$(9)

When $s_1=s_2=\dots =s_r=\mathit{p}$, the $\mathit{I}\mathit{R}{ }_r(\mathit{s},\mathrm{\Sigma })$ reduces to the inverse Wishart $\mathcal{I}{\mathcal{W}}_r(\mathit{p},\mathrm{\Sigma })$ given as following.

Definition 2.6.

An $\mathit{r}\times \mathit{r}$ positive definite symmetric matrix $\mathit{Y}$ is said to have a matrix-variate inverse Wishart distribution with shape parameter $\mathit{p}>\frac{\mathit{r}-1}{2}$ and scale parameter $\mathrm{\Sigma }$ in ${\mathrm{\Omega }}_r$, if its probability density function is specified by

(10) $\mathcal{I}{\mathcal{W}}_r(\mathit{p},\mathrm{\Sigma })=\frac{{(det\mathrm{\Sigma })}^{-\mathit{p}}}{2^{\mathit{rp}}{\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{p})}{e}^{-\frac{1}{2}⟨{\mathrm{\Sigma }}^{-1},Y^{-1}⟩}(det\mathit{Y}{)}^{-\mathit{p}-\frac{\mathit{r}+1}{2}}{\mathbf{1}}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{Y}).$(10)

In the following, we briefly review the simulation procedures for random Riesz and inverse Riesz matrices. According to Kessentini et al. (2020), ${\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })$ and $\mathcal{I}{\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })$ are generated using an extension of the Bartlett decomposition (Bartlett 1934). Let $\mathrm{T}$ be in ${\mathcal{T}}_r$ with entries $t_{\mathit{ij}}$ satisfying:

1. ${\mathrm{t}}_{\mathit{ij}}$ is standard normal for all $1\le \mathit{j}<\mathit{i}\le \mathit{r}$.

2. $({\mathrm{t}}_{\mathit{ii}}{)}_{1\le \mathit{i}\le \mathit{r}}$ is $\sqrt{\mathit{\gamma }(\frac{2s_i-\mathit{i}+1}{2},\frac{1}{2})}$, where $\mathit{\gamma }(\frac{2s_i-\mathit{i}+1}{2},\frac{1}{2})$ is the gamma distribution with shape parameter $\frac{2s_i-\mathit{i}+1}{2}$; $s_i>\frac{\mathit{i}-1}{2}$ for all $1\le \mathit{i}\le \mathit{r}$ and scale parameter equal to $\frac{1}{2}$.

3. The $({\mathrm{t}}_{\mathit{ij}}{)}_{1\le \mathit{j}\le \mathit{i}\le \mathit{r}}$‘s are independent.

The matrix $\mathrm{X}=\mathrm{T}{\mathrm{T}}^{\ast }$ has a Riesz ${\mathit{R}}_r(\mathit{s},I_r)$ distribution, where $\mathit{s}=(s_1,\dots ,s_r)$ and $I_r$ is the identity matrix. Handling a matrix other than the identity matrix is easy. If the $\mathit{r}\times \mathit{r}$ symmetric matrix $\mathrm{\Sigma }$ has the Cholesky factorization $\mathrm{\Sigma }=\mathit{L}{L}^{\ast }$, with $\mathit{L}\in {\mathcal{T}}_r$, then $\mathrm{Z}=\mathit{L}\mathrm{T}(\mathit{L}\mathrm{T}{)}^{\ast }$ has a ${\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })$. In fact, according to (Tounsi and Zine 2012), the matrix $\mathrm{Y}={\mathrm{Z}}^{-1}$ follows an $\mathcal{I}{\mathit{R}}_r(\mathit{s},\mathrm{\Sigma })$.

2.2. The mixture of the riesz probability distribution

In modern statistics, mixture models correspond to an important probabilistic model which are used for data classification providing efficient approaches for model-based clustering (Fraley and Raftery 2002). They have been developed so far by various authors like McLachlan et al. (2002), Melnykov and Maitra (2010) and Gormley et al. (2023). Hence, a wide mixture models examples can be reported including mixture of gamma, mixture of exponential, mixture of Gaussian and mixture of the Riesz distributions (Kessentini et al. 2020; He and Chen 2022; Shafiq et al. 2022; Huang and Gou 2023).

In the following subsection, we exhibit a description of the mixture of the Riesz probability distribution which corresponds to a weighted sum of $\mathit{K}\in 1\mathrm{N}$ component Riesz densities demonstrated given below.

Let ${\mathrm{X}}_1,{\mathrm{X}}_2,\dots ,{\mathrm{X}}_{\mathrm{N}}$ be $\mathit{N}$ i.i.d. matrices with common density mixture expressed by

(11) $\mathit{f}(X_i|\mathrm{\Phi })=\underset{\mathrm{k}=1}{\sum ^{\mathrm{K}}}{\pi }_{\mathrm{k}}{\mathrm{f}}_{\mathrm{k}}({\mathrm{X}}_{\mathrm{i}}|{\theta }_{\mathrm{k}}),$(11)

where ${\pi }_k$ represents the mixing weights, $0<{\pi }_k<1$ such that $\sum _{\mathit{k}=1}^K{\pi }_k=1$. The $f_k(.|{\theta }_k)$ denotes the Riesz probability density function with parameters ${\theta }_k=(s_k,{\mathrm{\Sigma }}_k)$, $\mathrm{\forall }1\le \mathit{k}\le \mathit{K}$ (see Equation (6)(6) ${\mathcal{R}}_r(\mathit{s},\mathrm{\Sigma })=\frac{1}{2^{\sum _{\mathit{j}=1}^r{s}_j}{\mathrm{\Gamma }}_{{\mathrm{\Omega }}_{\mathrm{r}}}(\mathit{s}){\mathrm{\Delta }}_s(\mathrm{\Sigma })}{e}^{-\frac{1}{2}⟨{\mathrm{\Sigma }}^{-1},\mathit{X}⟩}{\mathrm{\Delta }}_{\mathit{s}-\frac{\mathit{r}+1}{2}}(\mathit{X}){\mathbf{1}}_{{\mathrm{\Omega }}_r}(\mathit{X}),$(6) ), where $s_k=(s_{\mathit{k},1},s_{\mathit{k},2},\dots ,s_{\mathit{k},\mathit{r}})$. In addition, $\mathrm{\Phi }=({\pi }_1,\dots ,{\pi }_K,{\theta }_1,\dots {\theta }_K)$ denotes the vector of unknown parameters to be estimated by the EM algorithm. The standard reference on the EM algorithm and its convergence is Dempster et al. (1977).

At this stage of analysis, we briefly recall the principle of the use of the EM algorithm. It is often used to estimate parameters of a mixture model, in which the exact component model from which a data point is generated is hidden from us. It starts with randomly assigning values to all the parameters to be estimated. It then iterately alternates between two steps, called the Expectation step (E-step) and the Maximization step (M-step), respectively. In the E-step, it computes the expected likelihood for the complete data (the so-called $\mathbf{Q}$-function) where the expectation is taken for the computed conditional distribution of the latent variables given the current settings of parameters and our observed (incomplete) data. In the M-step, it re-estimates all the parameters by maximizing the $\mathbf{Q}$-function. Once we have a new generation of parameter values, we can repeat the E-step and another M-step. This process continues until the likelihood converges, i.e. reaching a local maxima. Intuitively, what EM does is to iteratively increase the data by guessing the values of the hidden variables and to re-estimate the parameters by assuming that the guessed values are the true values. Recently, Kessentini et al. (2020) developed an Expectation-Maximization algorithm to estimate the parameters of the Riesz mixture model. Therefore, the parameters estimation of ${\pi }_k,s_k$ and ${\mathrm{\Sigma }}_k$, $\mathrm{\forall }1\le \mathit{k}\le \mathit{K}$ are calculated, according to the authors Kessentini et al. (2020), using these points.

• The unclassified mixture $\mathrm{X}=({\mathrm{X}}_1,\dots ,{\mathrm{X}}_{\mathrm{N}})$ is augmented by the hidden $\mathit{N}-$dimentional vector $\mathrm{z}=({\mathrm{z}}_1,\dots ,{\mathrm{z}}_{\mathrm{N}})$ where ${\mathrm{z}}_{\mathrm{i}}=\mathit{k}$ means that ${\mathrm{X}}_{\mathrm{i}}$ is issued from the $k^{\mathit{th}}$ distribution. The complete-data likelihood and log likelihood are formulated as

  (12) ${\mathit{L}}_c(\mathrm{\Phi })=\underset{\mathrm{i}=1}{\prod ^{\mathrm{N}}}\underset{\mathrm{k}=1}{\sum ^{\mathrm{K}}}\mathbf{1}\underset{\{{\mathrm{z}}_{\mathrm{i}}=\mathrm{k}\}}{}{\pi }_{\mathrm{k}}{\mathrm{f}}_{\mathrm{k}}({\mathrm{X}}_{\mathrm{i}}|{\theta }_{\mathrm{k}}),$(12)
  (13) $log\left({\mathit{L}}_c(\mathrm{\Phi })\right)=\underset{\mathit{i}=1}{\sum ^N}\underset{\mathit{k}=1}{\sum ^K}\mathbf{1}\underset{\{{\mathrm{z}}_{\mathrm{i}}=\mathit{k}\}}{}log\left({\pi }_k{f}_k(X_i|{\theta }_k)\right).$(13)

• The expectation step is equivalent to estimate at iteration $(\mathit{t}+1)$ the conditional expectation $\mathbf{Q}\left(\mathrm{\Phi }|{\mathrm{\Phi }}^{(\mathrm{t})}\right)$ of the complete-data log likelihood given the previous value of parameters ${\mathrm{\Phi }}^{(\mathrm{t})}$

  (14) $\begin{array}{cc}\mathbf{Q}\left(\mathrm{\Phi }|{\mathrm{\Phi }}^{(\mathit{t})}\right)& ={\mathbf{E}}_{{\mathrm{\Phi }}^{(\mathrm{t})}}\left[log\left({\mathit{L}}_c(\mathrm{\Phi })\right)\right]\\ \mathit{\hspace{0.17em}}& =\underset{\mathit{i}=1}{\sum ^N}\underset{\mathit{k}=1}{\sum ^K}{\mathbf{E}}_{{\mathrm{\Phi }}^{(\mathrm{t})}}\left[z_i|X_i\right]log\left({\pi }_k^{(t)}{f}_k\left(X_i|{\theta }_k^{(t)}\right)\right)\\ \mathit{\hspace{0.17em}}& =\underset{\mathit{i}=1}{\sum ^N}\underset{\mathit{k}=1}{\sum ^K}{\tilde{p}}_{\mathit{ik}}log\left({\pi }_k^{(t)}{f}_k\left(X_i|{\theta }_k^{(t)}\right)\right),\end{array}$(14)
  where ${\tilde{p}}_{\mathit{ik}}$ is the posterior probability of the component membership $P_{{\mathrm{\Phi }}^{(\mathit{t})}}\left[{\mathrm{z}}_{\mathrm{i}}=\mathit{k}|X_i\right]$ at iteration $(\mathit{t}+1)$ and which is given by the Bayes rule
  (15) $\begin{array}{ccc}{\tilde{p}}_{\mathit{ik}}& =& \frac{{\mathit{P}}_{{\mathrm{\Phi }}^{(\mathit{t})}}\left[{\mathrm{z}}_{\mathrm{i}}=\mathit{k}\right]f_k\left(X_i|{\theta }_k^{(t)}\right)}{\mathit{f}(X_i|{\mathrm{\Phi }}^{(\mathit{t})})}\\ \mathit{\hspace{0.17em}}& =& \frac{{\pi }_k^{(t)}{f}_k\left(X_i|{\theta }_k^{(t)}\right)}{\underset{\mathit{l}=1}{\sum ^K}{\pi }_l^{(t)}{f}_l\left(X_i|{\theta }_l^{(t)}\right)}.\end{array}$(15)

• The maximisation step of the algorithm consists of estimating, at the $(\mathit{t}+1{)}^{\mathit{th}}$ iteration, the new $\hat{\mathrm{\Phi }}={\mathrm{\Phi }}^{(\mathrm{t}+1)}$ that maximises the expectation subject to the constraint $\underset{\mathit{k}=1}{\sum ^K}{\pi }_k=1$

  (16) $\hat{\mathrm{\Phi }}=arg\underset{\mathrm{\Phi }}{max}\left(\mathbf{Q}\left(\mathrm{\Phi }|{\mathrm{\Phi }}^{(\mathit{t})}\right)\right).$(16)
  Given the concavity of the log likelihood, solving the first order condition is sufficient.

Therefore, the maximisation step gives a new estimate of the unknown parameters $({\pi }_k,s_k,{\mathrm{\Sigma }}_k)$ for all $1\le \mathit{k}\le \mathit{K}$ as indicated below.

We have the new estimate of the mixing proportion ${\pi }_k$ is given by

(17) ${\hat{\pi }}_{\mathit{k}}=\frac{1}{N}\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}.$(17)

The new estimate of the components $s_{\mathit{k},\mathit{j}}$ of the vector $s_k$ is expressed by

(18) ${\hat{s}}_{\mathit{k},\mathit{j}}=\left\{\begin{array}{ccc}{\psi }^{-1}\left(-log(2{\mathrm{\Delta }}_1({\mathrm{\Sigma }}_k))+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left({\mathrm{\Delta }}_1(X_i)\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& \mathit{j}=1\\ \frac{\mathit{j}-1}{2}+{\psi }^{-1}\left(log\left(\frac{{\mathrm{\Delta }}_{\mathit{j}-1}({\mathrm{\Sigma }}_k)}{2{\mathrm{\Delta }}_j({\mathrm{\Sigma }}_k)}\right)+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left(\frac{{\mathrm{\Delta }}_j(X_i)}{{\mathrm{\Delta }}_{\mathit{j}-1}(X_i)}\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& 2\le \mathit{j}\le \mathit{r}\end{array}\right.,$(18)

where the utilised function $\mathit{\psi }$ is the digamma function defined as the logarithmic derivative of the gamma function. We have $\mathit{\psi }=\frac{{\mathrm{\Gamma }}^{\prime }}{\mathrm{\Gamma }}$.

And, finally the new estimate of the scale matrix ${\mathrm{\Sigma }}_k$ of the Riesz distribution is determined by

(19) $\underset{\mathit{j}=1}{\sum ^r}({\hat{s}}_{\mathit{k},\mathit{j}}-{\hat{s}}_{\mathit{k},\mathit{j}+1})\widehat{{\mathrm{\Sigma }}_k}{\left(P_j(\widehat{{\mathrm{\Sigma }}_k})\right)}^{-1}\widehat{{\mathrm{\Sigma }}_k}=\frac{1}{2\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}{X}_i,$(19)

where ${\hat{s}}_{\mathit{k},\mathit{r}+1}=s_{\mathit{k},\mathit{r}+1}=0.$

Note that the latter equation leads to a system of linear equations that can be solved, as described on the work done by Kessentini et al. (2020), to determine $\widehat{{\mathrm{\Sigma }}_k}$.

The EM algorithm for the Riesz mixture model is then summarised by Algorithm 1.

Table

Algorithm 1. The EM algorithm for the Riesz mixture model
1: begin
2: Initialisation: ${\mathrm{\Phi }}^{(0)}=\left({{\pi }_k}^{(0)},{s_k}^{(0)},{{\mathrm{\Sigma }}_k}^{(0)}\right)$
3: Iteration (t+$1):$
4: Expectation:
5: ${\tilde{p}}_{\mathit{ik}}=\frac{{\pi }_k^{(t)}{f}_k\left(X_i|{s_k}^{(0)},{{\mathrm{\Sigma }}_k}^{(0)}\right)}{\underset{\mathit{j}=1}{\sum ^K}{\pi }_j^{(t)}{f}_j\left(X_i|{s_j}^{(0)},{{\mathrm{\Sigma }}_j}^{(0)}\right)}.$
6: Maximization:
7: ${\hat{\pi }}_{\mathit{k}}=\frac{1}{N}\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}.$
8: ${\hat{s}}_{\mathit{k},\mathit{j}}=\left\{\begin{array}{ccc}{\psi }^{-1}\left(-log(2{\mathrm{\Delta }}_1({\mathrm{\Sigma }}_k))+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left({\mathrm{\Delta }}_1(X_i)\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& \mathit{j}=1\\ \frac{\mathit{j}-1}{2}+{\psi }^{-1}\left(log\left(\frac{{\mathrm{\Delta }}_{\mathit{j}-1}({\mathrm{\Sigma }}_k)}{2{\mathrm{\Delta }}_j({\mathrm{\Sigma }}_k)}\right)+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left(\frac{{\mathrm{\Delta }}_j(X_i)}{{\mathrm{\Delta }}_{\mathit{j}-1}(X_i)}\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& 2\le \mathit{j}\le \mathit{r}\end{array}\right.$
9: $\underset{\mathit{l}=1}{\sum ^r}({\hat{s}}_{\mathit{k},\mathit{l}}-{\hat{s}}_{\mathit{k},\mathit{l}+1})\widehat{{\mathrm{\Sigma }}_k}{\left(P_l(\widehat{{\mathrm{\Sigma }}_k})\right)}^{-1}\widehat{{\mathrm{\Sigma }}_k}=\frac{1}{2\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}{X}_i.$
10: If $||{\mathrm{\Phi }}^{(\mathit{t}+1)}-{\mathrm{\Phi }}^{(\mathit{t})}||<\mathit{\epsilon }$ is not satisfied return to step $2.$
11: end

3. Application

This section is devoted to the application of the mixture of Riesz probability distributions for medical image segmentation. The significance of the latter refers to its ability to extract information about the human body, and more specifically about tissue characterisation, organs, anatomic structures, lesions and tumours. This information helps doctors delineate and track the progress of diseases. The brain is a particulary complex structure, and its segmentation is an important step for solving many problems. We are basically interested in brain MR image unsupervised segmentation. Our objective is therefore, to segment brain tumours by using the proposed Riesz mixture model for both simulated image and real brain MR image in order to identify two regions (tumour and non-tumour). The creation of a robust and a flexible algorithm to ensure the best image segmentation result based on Riesz distribution is hence necessary. In multiple classical works, the segmentation process rests on the use of the Gaussian distribution (Fenyi et al. 2020; Yang et al. 2022). However, in many fields like finance, handwriting recognition or image segmentation, mixture models were commonly extended from Gaussian mixture model to the multivariate one in order to fit vectors of unknown parameters. They were even extended to fit matrix samples as shown by Hidot and Saint-Jean using the Wishart probability distribution (Hidot and Saint-Jean 2010).

Until achieving the final segmentation result using Riesz mixture model and verifying the high performance of this algorithm over than the Wishart one, such stages must be followed (see Figure 1 below).

Figure 1. The flowchart of the proposed image segmentation method.

The proposed algorithm of segmentation relies upon eight steps, which are described as follows.

• Step 1: Divided image into blocks and choice of the number of image regions $(\mathit{K}=2)$ for each block.

• Step 2: Input observed image in a vector of the matrix of pixels.

• Step 3: Initialization of the vector of the parameters $\mathrm{\Phi }=\left({\pi }_k,s_k,{\mathrm{\Sigma }}_k\right)$, $\mathrm{\forall }\mathit{k}=1,2$ by using the Fuzzy C-means algorithm (Dunn 1973; Yu et al. 2023).

• Step 4 (Expectation): Computing the posterior probabilities at iteration $(\mathit{t}+1)$ by the Bayes rule

  ${\tilde{p}}_{\mathit{ik}}=\frac{{\pi }_k^{(t)}{f}_k\left(X_i|{s_k}^{(0)},{{\mathrm{\Sigma }}_k}^{(0)}\right)}{\underset{\mathit{j}=1}{\sum ^K}{\pi }_j^{(t)}{f}_j\left(X_i|{s_j}^{(0)},{{\mathrm{\Sigma }}_j}^{(0)}\right)},\mathrm{for}\mathrm{all}\mathit{k}=1,2.$

• Step 5 (Maximization): Maximizing the conditional expectation

  $\mathbf{Q}\left(\mathrm{\Phi }|{\mathrm{\Phi }}^{(\mathit{t})}\right)=\underset{\mathit{i}=1}{\sum ^N}\underset{\mathit{k}=1}{\sum ^K}{\tilde{p}}_{\mathit{ik}}log\left({\pi }_k^{(t)}{f}_k\left(X_i|s_k^{(t)},{\mathrm{\Sigma }}_k^{(t)}\right)\right),$
  with respect to $s_k,{\mathrm{\Sigma }}_k$ and ${\pi }_k$, $\mathrm{\forall k}=1,2$. The update of the Riesz mixture parameters (the mixing weight ${\pi }_k\in \left]0,1\right[$, the shape parameter of the Riesz distribution $s_k=(s_{\mathit{k},1},s_{\mathit{k},2},\dots ,s_{\mathit{k},\mathit{r}})\in {\mathbb{R}}^r$, where $s_{\mathit{k},\mathit{j}}>\frac{\mathit{j}-1}{2},\mathrm{\forall }1\le \mathit{j}\le \mathit{r}$ and its scale parameter ${\mathrm{\Sigma }}_k\in {\mathrm{\Omega }}_r$) is supplied by
  ${\hat{\pi }}_{\mathit{k}}=\frac{1}{N}\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}},$
  ${\hat{s}}_{\mathit{k},\mathit{j}}=\left\{\begin{array}{ccc}{\psi }^{-1}\left(-log(2{\mathrm{\Delta }}_1({\mathrm{\Sigma }}_k))+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left({\mathrm{\Delta }}_1(X_i)\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& \mathit{j}=1\\ \frac{\mathit{j}-1}{2}+{\psi }^{-1}\left(log\left(\frac{{\mathrm{\Delta }}_{\mathit{j}-1}({\mathrm{\Sigma }}_k)}{2{\mathrm{\Delta }}_j({\mathrm{\Sigma }}_k)}\right)+\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}log\left(\frac{{\mathrm{\Delta }}_j(X_i)}{{\mathrm{\Delta }}_{\mathit{j}-1}(X_i)}\right)}{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)& ;& 2\le \mathit{j}\le \mathit{r}\end{array}\right.,$
  In the case $\mathit{r}=3,$ the rows of the symmetric matrix ${\hat{\mathrm{\Sigma }}}_{\mathit{k}}$ can be computed in the order following
  ${\hat{\mathrm{\Sigma }}}_{k_{1\mathit{j}}}=\left(\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}{X}_{{\mathit{i}}_{1\mathit{j}}}}{2\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}\right)/{\hat{s}}_{\mathit{k},1},\mathrm{for}\mathit{j}=1,2,3.$
  then
  (20) ${\hat{\mathrm{\Sigma }}}_{k_{2\mathit{j}}}=\left(\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}{X}_{{\mathit{i}}_{2\mathit{j}}}}{2\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}-({\hat{s}}_{\mathit{k},1}-{\hat{s}}_{\mathit{k},2})\frac{{\hat{\mathrm{\Sigma }}}_{k_{12}}{\hat{\mathrm{\Sigma }}}_{k_{1\mathit{j}}}}{{\hat{\mathrm{\Sigma }}}_{k_{11}}}\right)/{\hat{s}}_{\mathit{k},2},\mathrm{for}\mathit{j}=2,3.$(20)
  and finally
  (21) ${\hat{\mathrm{\Sigma }}}_{k_{33}}=\left(\frac{\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}{X}_{{\mathit{i}}_{33}}}{2\underset{\mathit{i}=1}{\sum ^N}{\tilde{p}}_{\mathit{ik}}}-({\hat{s}}_{\mathit{k},1}-{\hat{s}}_{\mathit{k},2})\frac{{({\hat{\mathrm{\Sigma }}}_{k_{13}})}^2}{{\hat{\mathrm{\Sigma }}}_{k_{11}}}-({\hat{s}}_{\mathit{k},2}-{\hat{s}}_{\mathit{k},3})\frac{{\hat{\mathrm{\Sigma }}}_{k_{22}}{({\hat{\mathrm{\Sigma }}}_{k_{13}})}^2-2{\hat{\mathrm{\Sigma }}}_{k_{12}}{\hat{\mathrm{\Sigma }}}_{k_{13}}{\hat{\mathrm{\Sigma }}}_{k_{23}}+{\hat{\mathrm{\Sigma }}}_{k_{11}}{({\hat{\mathrm{\Sigma }}}_{k_{23}})}^2}{{\hat{\mathrm{\Sigma }}}_{k_{11}}{\hat{\mathrm{\Sigma }}}_{k_{22}}-{\hat{\mathrm{\Sigma }}}_{k_{12}}^2}\right)/{\hat{s}}_{\mathit{k},3}.\mathit{a}$(21)
  More details are given by Kessentini et al. (2020).

• Step 6: Iterate steps 4 and 5 until the following constraint is obtained

  (22) $||{\mathrm{\Phi }}^{(\mathit{t}+1)}-{\mathrm{\Phi }}^{(\mathit{t})}||<\mathit{ϵ}.$(22)
  That is the convergence of the EM algorithm for the Riesz mixture model (EMR) is assumed when the absolute difference between successive estimates is less than ${10}^{-3}$. Therefore, the stop criterion of the EMR is accomplished after 20 iterations.

• Step 7: Classify all the matrix of pixels into one of the $\mathit{K}$ classes by calculating the posterior probability denoted $\mathbb{P}(\mathit{\alpha };s_k,{\mathrm{\Sigma }}_k)$ of the matrix of pixels $\mathit{\alpha }$ for each class $C_1$ and $C_2$. The image pixels are then classified (labelled) based on the highest posterior probability. Hence, the matrix of pixels denoted $\mathit{\alpha }$ is said to be a tumour if

  (23) $\mathbb{P}(\mathit{\alpha };s_1,{\mathrm{\Sigma }}_1)>\mathbb{P}(\mathit{\alpha };s_2,{\mathrm{\Sigma }}_2);$(23)
  otherwise, it is considered as normal tissue.

• Step 8: Construct labelled image corresponding of each block.

In order to prove further the effectiveness of the proposed algorithm and examine the quality of our segmentation results, some performance metrics are defined in the following subsection.

3.1. Performance metrics

Resting on the original image $\mathit{I}$ and the segmented image $\hat{I}$, six evaluation criteria can be calculated: $\mathit{PSNR}$, $\mathit{MCR}$, $\mathit{RI}$, $\mathit{ARI}$, $\mathit{MI}$ and $\mathit{HI}$. These performance measure parameters can be found in (Saladi and Amutha Prabha 2017; Fouad et al. 2017; Correa-Morris et al. 2022; Zitouni et al. 2023; Sundqvist et al. 2023) and (Somashekara and Manjunatha 2014).

The Peak Signal to Noise Ratio $(\mathit{PSNR})$ in $\mathit{dB}$ is most commonly used as a quality measure of reconstruction that is expressed by

(24) $\mathit{PSNR}=10log\left(\frac{Ima{x}^2}{\mathit{MSE}}\right),$(24)

where, $\mathit{Imax}$ corresponds to the maximum intensity of pixel value of the image and $\mathit{MSE}$ represents the Mean Squared Error. Small $\mathit{MSE}$ value means less errors, resulting in high $\mathit{PSNR}$. Acceptable $\mathit{PSNR}$ range is $>20$ $\mathit{dB}$.

In order to measure the segmentation accuracy, we also need to define the Misclassification Ratio $(\mathit{MCR})$

(25) $\mathit{MCR}=\frac{\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{misclassified}\mathrm{pixels}}{\mathrm{the}\mathrm{number}\mathrm{of}\mathrm{all}\mathrm{pixels}\mathrm{from}\mathrm{the}\mathrm{original}\mathrm{image}}.$(25)

We have $\mathit{MCR}\in [0,1]$. The smaller $\mathit{MCR}$ is the better the segmentation result becomes.

To further assess the performances of the proposed algorithm, we used other image segmentation evaluation metrics which are indicated as follows.

Rand Index $(\mathit{RI})$ is a measure of the similarity between two data clusterings. It has a value between 0 and 1, with 0 indicating that the two data clusterings do not agree on any pair of points and 1 indicating that the data clusterings are exactly the same.

The Adjusted Rand Index $(\mathit{ARI})$ is defined as the similarity between the actual labels and predicted labels, ignoring permutations. The $\mathit{ARI}$ has a range of $[-1,1],$ negative values are considered as bad clusterings and the value one means perfect clustering. The $\mathit{ARI}$ is denoted as

(26) $\mathit{ARI}=\frac{\mathit{RI}-\mathit{E}[\mathit{RI}]}{\mathit{max}(\mathit{RI})-\mathit{E}[\mathit{RI}]},$(26)

where $\mathit{E}[\mathit{RI}]$ is the expected value of $\mathit{RI}$.

Further statistical tests for clustering evaluation are also depicted below.

The Mirkin Metric $(\mathit{MI})$ is a transformation of Rand Index. It is the probability of a disagreement expressed by

(27) $\mathit{MI}=\mathit{n}(\mathit{n}-1)(\mathit{RI}-1),$(27)

where $\mathit{n}$ is the number of data items.

The Hubert index $(\mathit{HI})$ is the difference between the probability of an agreement and a disagreement.

These measures are helpful to evaluate the ability of the proposed segmentation algorithm and discriminate tissues especially tumours.

3.2. Experimental results

This section presents the experimental results of the proposed segmentation algorithm with synthetic and real brain MR image.

3.2.1. Synthetic image segmentation

At this stage of analysis, we shall address the obtained results, discuss them and evaluate them in order to assess the effectiveness and feasibility of our approach. The proposed method is evaluated first using synthetic image. A simulated image is generated by the use of the Riesz distributions. Using Matlab software, simulated image $\hat{I}$ is generated by a sample of total size $8000$ of random matrices. Each matrix belongs to the cone of $3\times 3$ symmetric positive definite matrices denoted by ${\mathrm{\Omega }}_3$. Thus, we generate the pixels of these matrices based on the extension of the Bartlett decomposition for the Riesz distribution (as described in Section 2) with gamma on the diagonal and Gaussian off-diagonal. The simulated image $\hat{I}$ is composed of $\mathit{K}=2$ classes as described below.

We consider a mixture of two Riesz probability distributions (Riesz mixture model 1) generated using the extended Bartlett decomposition of Riesz random matrices (Kessentini et al. 2020) on ${\mathrm{\Omega }}_3$. From this perspective, let

• ${\mathrm{M}}_1\sim {\mathcal{R}}_3(s_1,{\mathrm{\Sigma }}_1)$, where $s_1=(2,5,6)$ and ${\mathrm{\Sigma }}_1=\left(\begin{array}{ccc}2.2& 0.2& 0.2\\ 0.2& 2.2& 0.2\\ 0.2& 0.2& 2.2\end{array}\right).$

• ${\mathrm{M}}_2\sim {\mathcal{R}}_3(s_2,{\mathrm{\Sigma }}_2)$, where $s_2=(3,4,7)$ and ${\mathrm{\Sigma }}_2=\left(\begin{array}{ccc}4.3& 0.3& 0.3\\ 0.3& 4.3& 0.3\\ 0.3& 0.3& 4.3\end{array}\right).$

Recall that, to generate the two random matrices ${\mathrm{M}}_1$ and ${\mathrm{M}}_2$, we set $\mathrm{X}={\mathrm{T}}_1{\mathrm{T}}_1^{\ast }$ and $\mathrm{Y}={\mathrm{T}}_2{\mathrm{T}}_2^{\ast }$, where ${\mathrm{T}}_1$ and ${\mathrm{T}}_2$ are in ${\mathcal{T}}_r$ with independent elements such that ${\mathrm{T}}_k=({\mathrm{t}}_{{\mathit{k}}_{\mathit{ij}}}{)}_{1\le \mathit{j}\le \mathit{i}\le \mathit{r}}$ for all $\mathit{k}=1,2$. Thus, we generate the diagonal elements ${\mathrm{t}}_{{\mathit{k}}_{\mathit{ii}}}^2$ with the usual real gamma distribution $\mathit{\gamma }(\frac{2{\mathit{s}}_{\mathit{k},\mathit{i}}-\mathit{i}+1}{2},\frac{1}{2})$; $s_{\mathit{k},\mathit{i}}>\frac{\mathit{i}-1}{2}$ for all $1\le \mathit{i}\le \mathit{r}$ and we generate the off-diagonal elements ${\mathrm{t}}_{{\mathit{k}}_{\mathit{ij}}}$ with the Gaussian distribution $\mathcal{N}(0,1)$; $\mathrm{\forall }1\le \mathit{j}<\mathit{i}\le \mathit{r}$. Therefore, we get $\mathrm{X}\sim R_3(s_1,I_3)$ and $\mathrm{Y}\sim R_3(s_2,I_3)$. Now, by using the Cholesky decomposition of ${\mathrm{\Sigma }}_k=L_k{L}_k^{\ast }$, where $L_k\in {\mathcal{T}}_r$ for all $\mathit{k}=1,2$, we construct the random matrices ${\mathrm{M}}_k=L_k{\mathrm{T}}_k{\left(L_k{\mathrm{T}}_k\right)}^{\ast }$ which follow the Riesz distribution $R_3(s_k,{\mathrm{\Sigma }}_k)$.

Besides, for further statistical analysis and simulation studies, we introduce other mixture of Riesz probability distributions with other parameters which are given in Table 1 below. Note that we adopt the same scale matrices as in Ref (Hidot and Saint-Jean 2010). for the third Riesz mixture model.

Table 1. Other riesz mixture models.

Riesz mixture model 2	${\mathrm{M}}_1^{\prime }\sim {\mathcal{R}}_3(s_1^{\prime },{\mathrm{\Sigma }}_1^{\prime })$ and ${\mathrm{M}}_2^{\prime }\sim {\mathcal{R}}_3(s_2^{\prime },{\mathrm{\Sigma }}_2^{\prime })$
$s_1^{\prime }$	$(4.4151,2.1242,6.7157)$
$s_2^{\prime }$	$(4.3197,2.1438,7.3523)$
${\mathrm{\Sigma }}_1^{\prime }$	$\left(\begin{array}{ccc}1.1346& 0.0658& 0.0519\\ 0.0658& 1.3122& -0.0304\\ 0.0519& -0.0304& 2.5979\end{array}\right)$
${\mathrm{\Sigma }}_2^{\prime }$	$\left(\begin{array}{ccc}1.1606& 0.0607& 0.0771\\ 0.0607& 1.2805& 0.1208\\ 0.0771& 0.1208& 2.3083\end{array}\right)$
Riesz mixture model 3	${\mathrm{M}}_1^{{ }^{\prime \prime }}\sim {\mathcal{R}}_3(s_1^{{ }^{\prime \prime }},{\mathrm{\Sigma }}_1^{{ }^{\prime \prime }})$ and ${\mathrm{M}}_2^{{ }^{\prime \prime }}\sim {\mathcal{R}}_3(s_2^{{ }^{\prime \prime }},{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }})$
$s_1^{{ }^{\prime \prime }}$	$(10,17,19)$
$s_2^{{ }^{\prime \prime }}$	$(15,17,18)$
${\mathrm{\Sigma }}_1^{{ }^{\prime \prime }}$	$\left(\begin{array}{ccc}2.8255& 1.6921& -2.9052\\ 1.6921& 2.6183& -2.2620\\ -2.9052& -2.2620& 3.7281\end{array}\right)$
${\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}$	$\left(\begin{array}{ccc}2.9232& 1.8945& -0.4002\\ 1.8945& 2.7291& 0.7391\\ -0.4002& 0.7391& 1.4830\end{array}\right).$

To this extent, we generate a sample ${\mathrm{X}}_1,\dots ,{\mathrm{X}}_{\mathrm{N}}$ of total size $\mathit{N}=8000$ belonging to one of the two Riesz probability distributions for each Riesz mixture model, with mixing weights for the first class $C_1$ and the second class $C_2$ that are indicated respectively as ${\pi }_1=0.4$ and ${\pi }_2=0.6$. To check whether we have generated matrix of Riesz distribution, comparative studies between the empirical mean and the expected value of the Riesz random matrix (see Equation (8)(8) $\mathit{E}(\mathit{X})=-\sum _{\mathit{l}=1}^r(s_{\mathit{r}-\mathit{l}}-s_{\mathit{r}-\mathit{l}+1}){\left(P_l^{\ast }(\mathrm{\Sigma })\right)}^{-1}.$(8) )are carried out. The parameter estimates are obtained by the EM algorithm and the convergence is assumed when the absolute differences between successive estimates are less than ${10}^{-3}.$ So, after $20$ iterations, the task is accomplished.

In order to prove the efficiency and the superiority of the EM algorithm for the Riesz mixture models (EMR) over the EM algorithm for the Wishart mixture models (EMW), and examine the quality of our estimation results, we can adopt the following two metrics.

The Euclidean distance $\mathit{d}(s_k,\widehat{s_k})=||s_k-\widehat{s_k}|{|}_2$ for the shape parameter $s_k$, and the Hilbert-Schmidt distance given by

(28) $d_{\mathit{HS}}({\mathrm{\Sigma }}_k,\widehat{{\mathrm{\Sigma }}_k})=\sqrt{\mathrm{tr}\left(({\mathrm{\Sigma }}_k-\widehat{{\mathrm{\Sigma }}_k})({\mathrm{\Sigma }}_k-\widehat{{\mathrm{\Sigma }}_k}{)}^{\ast }\right)},$(28)

for the scale matrix ${\mathrm{\Sigma }}_k$.

The shape parameters and the scale matrices estimates for the different Riesz mixture models and their distances are plotted in Table 2 below.

Table 2. The resulting estimated parameters using the EMR and their distances toward the actual parameters.

Parameter	initial parameter value	Estimated value	Distance
$\widehat{s_1}$	$(1.5,3.2,4)$	$(1.8023,3.9629,4.5907)$	$\mathit{d}(s_1,\widehat{s_1})=1.7609$
$\widehat{s_2}$	$(2.3,2.9,4.5)$	$(2.7748,3.3752,4.9554)$	$\mathit{d}(s_2,\widehat{s_2})=2.1117$
$\widehat{{\mathrm{\Sigma }}_1}$	$\left(\begin{array}{ccc}1.5& 0.1& 0.11\\ 0.1& 1.8& 0.15\\ 0.11& 0.15& 2\end{array}\right)$	$\left(\begin{array}{ccc}1.9737& 0.1515& 0.1356\\ 0.1515& 2.0178& 0.1723\\ 0.1356& 0.1723& 2.1898\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_1,\widehat{{\mathrm{\Sigma }}_1})=0.3147$
$\widehat{{\mathrm{\Sigma }}_2}$	$\left(\begin{array}{ccc}2.5& 0.2& 0.21\\ 0.2& 3& 0.23\\ 0.21& 0.23& 4\end{array}\right)$	$\left(\begin{array}{ccc}2.9821& 0.2478& 0.2579\\ 0.2478& 3.9992& 0.2555\\ 0.2579& 0.2555& 4.2887\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2,\widehat{{\mathrm{\Sigma }}_2})=1.3566$
$\widehat{s_1^{\prime }}$	$\left(4.2,2,6.5\right)$	$\left(4.6827,2.1274,6.8422\right)$	$\mathit{d}(s_1^{\prime },\widehat{s_1^{\prime }})=0.2960$
$\widehat{s_2^{\prime }}$	$\left(4,2,7\right)$	$\left(4.5869,2.1513,7.5185\right)$	$\mathit{d}(s_2^{\prime },\widehat{s_2^{\prime }})=0.3147$
$\widehat{{\mathrm{\Sigma }}_1^{\prime }}$	$\left(\begin{array}{ccc}1& 0.08& 0.1\\ 0.08& 1.5& 0.1\\ 0.1& 0.1& 2.1\end{array}\right)$	$\left(\begin{array}{ccc}1.0686& 0.0633& 0.0297\\ 0.0633& 1.2998& -0.0600\\ 0.0297& -0.0600& 2.5521\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_1^{\prime },\widehat{{\mathrm{\Sigma }}_1^{\prime }})=0.0967$
$\widehat{{\mathrm{\Sigma }}_2^{\prime }}$	$\left(\begin{array}{ccc}1& 0.02& 0.07\\ 0.02& 1.5& 0.5\\ 0.07& 0.5& 2\end{array}\right)$	$\left(\begin{array}{ccc}1.0910& 0.0557& 0.0591\\ 0.0557& 1.2638& 0.1191\\ 0.0591& 0.1191& 2.2518\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2^{\prime },\widehat{{\mathrm{\Sigma }}_2^{\prime }})=0.0950$
$\widehat{s_1^{{ }^{\prime \prime }}}$	$\left(9,15,18.5\right)$	$\left(9.7831,16.8324,19.1857\right)$	$\mathit{d}(s_1^{{ }^{\prime \prime }},\widehat{s_1^{{ }^{\prime \prime }}})=0.3311$
$\widehat{s_2^{{ }^{\prime \prime }}}$	$\left(13,16.5,19\right)$	$\left(14.0143,16.9423,19.0237\right)$	$\mathit{d}(s_2^{{ }^{\prime \prime }},\widehat{s_2^{{ }^{\prime \prime }}})=1.4223$
$\widehat{{\mathrm{\Sigma }}_1^{{ }^{\prime \prime }}}$	$\left(\begin{array}{ccc}2.3& 1.5& 0.5\\ 1.5& 2& -2\\ 0.5& -2& 3\end{array}\right)$	$\left(\begin{array}{ccc}2.6961& 1.8729& -2.9499\\ 1.8729& 2.8549& -2.6202\\ -2.9499& -2.6202& 3.5011\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_1^{{ }^{\prime \prime }},\widehat{{\mathrm{\Sigma }}_1^{{ }^{\prime \prime }}})=0.6710$.
$\widehat{{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}}$	$\left(\begin{array}{ccc}2& 1.5& -0.2\\ 1.5& 3& 0.7\\ -0.2& 0.7& 1.8\end{array}\right)$	$\left(\begin{array}{ccc}2.5720& 1.9035& -0.5124\\ 1.9035& 2.9114& 0.5881\\ -0.5124& 0.5881& 1.6921\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2^{{ }^{\prime \prime }},\widehat{{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}})=0.5208$.

In the following (Table 3), we shall compare the EMR with the EMW.

Table 3. The resulting estimated parameters using the EMW and their distances toward the actual parameters.

Parameter	Estimated value	Distance
$\widehat{{\mathrm{\Sigma }}_1}$	$\left(\begin{array}{ccc}1.9097& 1.0476& 1.2616\\ 1.0476& 2.3129& 0.7130\\ 1.2616& 0.7130& 2.4455\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_1,\widehat{{\mathrm{\Sigma }}_1})=2.1165$
$\widehat{{\mathrm{\Sigma }}_2}$	$\left(\begin{array}{ccc}2.9905& 0.2688& 0.7079\\ 0.2688& 2.8791& 0.8410\\ 0.7079& 0.8410& 3.0745\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2,\widehat{{\mathrm{\Sigma }}_2})=2.4811$
$\widehat{{\mathrm{\Sigma }}_1^{\prime }}$	$\left(\begin{array}{ccc}1.0696& 0.0563& 0.0487\\ 0.0563& 1.3157& 0.0193\\ 0.0487& 0.0193& 2.5297\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_1^{\prime },\widehat{{\mathrm{\Sigma }}_1^{\prime }})=0.1184$
$\widehat{{\mathrm{\Sigma }}_2^{\prime }}$	$\left(\begin{array}{ccc}1.0934& 0.0507& 0.0736\\ 0.0507& 1.2819& 0.1843\\ 0.0736& 0.1843& 2.2260\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2^{\prime },\widehat{{\mathrm{\Sigma }}_2^{\prime }})=0.1399$
$\widehat{{\mathrm{\Sigma }}_1^{{ }^{\prime \prime }}}$	$\left(\begin{array}{ccc}2.8907& 2.0150& -1.9834\\ 2.0150& 2.9114& -1.7598\\ -1.9834& -1.7598& 3.8902\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2^{{ }^{\prime \prime }},\widehat{{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}})=1.5902$
$\widehat{{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}}$	$\left(\begin{array}{ccc}2.9971& 1.0032& -0.4358\\ 1.0032& 3.0415& 0.8120\\ -0.4358& 0.8120& 3.8785\end{array}\right)$	$d_{\mathit{HS}}({\mathrm{\Sigma }}_2^{\prime {\mathit{ }}^{\mathrm{\prime }}},\widehat{{\mathrm{\Sigma }}_2^{{ }^{\prime \prime }}})=1.4223$.

Referring to Tables 2 and 3, we assume that the lower Euclidean distance and $d_{\mathit{HS}}$ are reached by the EMR algorithm. Hence, it is inferred that the EMR performs better than the EMW. This shows that the Riesz probability distributions would even have more applications than the Wishart ones as they have more parameters (the shape parameter has a distant components) and then offer a higher degree of freedom. As we described above, simulation is used to generate synthetic images (See Figure 2 below).

Figure 2. The original synthetic image.

Therefore, in order to validate the proposed EMR algorithm, the latter was applied with synthetic images. The synthetic images associated to the three mixture Riesz models consists of two strongly noisy regions ($\mathit{K}=2$). It is portrayed in Figure 3.

Figure 3. Synthetic noisy images for the different riesz mixture models.

Various kinds of statistical methods have been elaborated for image segmentation. Among them, we mention those that apply Markov random field (MRF) model (Mouna et al. 2020). For instance, the region-based classification of polarimetric Synthetic Aperture Radar (SAR) images which uses Wishart MRF (Yin et al. 2020). The segmentation strategy based on the Wishart distribution is implemented for comparison. Owing to the best form of the Riesz distribution, more reliable results would be obtained. The corresponding segmentation results are illustrated in Figures 4–6 below.

Figure 4. Simulated segmented images by the mixture model 1.

Figure 5. Simulated segmented images by the mixture model 2.

Figure 6. Simulated segmented images by the mixture model 3.

Although there is much noise, we may record a great possibility of identifying the class of each pixel using the proposed Riesz probability distribution rather than the classical method based on Wishart probability distribution. Therefore, good image segmentation results are obtained using the new strategy through the first step estimation of the parameters. The estimators are very close to the true parameters values. Note that the estimated shape parameter by EMW is almost the average value of the associated vector $s_k$. Referring to the segmented images (Figures 4–6), the Riesz distribution can better characterise the class in heterogeneous areas and distinguish the regions with different textural information.

The table presented below (Table 4) outlines the quantitative results of the synthetic images.

Table 4. Different evaluation criteria values with Riesz and Wishart distributions.

Mixture model	$\mathit{MCR}$	$\mathit{PSNR}$	$\mathit{ARI}$	$\mathit{RI}$	$\mathit{MI}$	$\mathit{HI}$
Mixture Riesz model 1	$0.0035$	$24.3128$	$0.0610$	$0.6927$	$0.2587$	$0.3767$
Mixture Wishart model 1	$0.0052$	$22.1668$	$0.0579$	$0.6884$	$0.2642$	$0.3854$
Mixture Riesz model 2	$0.0198$	$21.0030$	$0.7049$	$0.8577$	$0.1423$	$-0.4183$
Mixture Wishart model 2	$0.1132$	$18.7668$	$0.0574$	$0.3648$	$0.6352$	$-0.2703$
Mixture Riesz model 3	$0.0455$	$21.6783$	$0.7028$	$0.8986$	$0.1014$	$0.2139$
Mixture Wishart model 3	$0.0628$	$20.9922$	$0.2471$	$0.6070$	$0.3930$	$0.4021$

Table 4 helps us derive the average values of MCR, $\mathit{PSNR},\mathit{ARI},\mathit{RI},\mathit{MI}$ and $\mathit{HI}$. The lowest $\mathit{MCR}$ and the highest $\mathit{PSNR}$ values illustrate the good performance of the proposed clustering method, which implies less errors and higher quality of the segmented image. The $\mathit{ARI}$ and the $\mathit{RI}$ values obtained of the Riesz distributions are higher than those obtained with the Wishart distributions. The lowest $\mathit{MI}$ and $\mathit{HI}$ values also indicate that the proposed algorithm based on the Riesz is more suitable than the Wishart one, which corroborates the effectiveness as well as reliability of the proposed algorithm. Between Riesz and Wishart distributions, there exists a tiny variation but as the different metrics values are more accepted with the new strategy, we can consider that the latter is much better to classify tissues than the classical method. The $\mathit{RI}$ values are, respectively, $0.6927$, $0.8577$ and $0.8986$ using the Riesz distributions. There are close to $1$ and are higher than the ones obtained by the Wishart distributions. From this perspective, the most correct and plausible segmentation result is provided by the proposed technique. The experiments reveal that these settings work well and can achieve satisfactory results.

3.2.2. Brain image segmentation

In this subsection, we highlight and we further discuss the robustness of the proposed method for settling brain image segmentation problem.

Firstly, we assume that the pixels intensities are independent samples from a mixture of the Riesz probability distributions. In order to illustrate the proposed distribution and the algorithm of segmentation, we consider the original real brain MR image and its blocks (see Figure 7). The image is of size ($183\times 183$) given from the Magnetic Resonance Imaging (MRI) Center of Sfax, Road of El-An, Sfax-Tunisia. In order to segment the brain image, we have divided the image into blocks.

Figure 7. Real brain MR image and its blocks.

We modelled each block of the brain image by a Riesz mixture model. We applied the proposed method to each block. Therefore, we suppose that each block of real brain MR image can be divided into two classes $(\mathit{K}=2)$: class of tumor and class of normal tissue. Each block is associated with mixing weights ${\pi }_k$, which should satisfy ${\sum }_{k=1}^2{\pi }_k=1$. To classify the matrix of pixels, we shall first estimate the mixture model parameters by using EM algorithm where initialisation of parameters are satisfied using the Fuzzy C-means algorithm (Kannan et al. 2010; Chaira 2011). Secondly, the classifier labels a test sample with the class that has the highest likelihood value. The classification result of pixels matrix for each block proved to be accurate and yielded satisfactory results.

In order to visualise the difference between the proposed segmentation algorithm with Riesz distribution and the proposed segmentation algorithm with Wishart distribution, the proposed technique has been applied on various brain images volume (See Figures 8 and 9).

Figure 8. Segmentation results of the brain image. (a) the segmented image by the riesz distribution, (b) the segmented image by the Wishart distribution.

Figure 9. Segmentation results of the brain image. (c) the segmented image by the riesz distribution, (d) the segmented image by the Wishart distribution.

The segmentation results are plotted in Figures 8 and 9 using the Riesz and the Wishart distributions, respectively. It is worth noting that the selection of the exact tumour region is intrinsic for a successful surgical intervention.

Referring to the segmented images by the Riesz (Figures 8(a) and 9(c) with respect to the segmented image by the Wishart (Figures 8(b) and 9(d), we infer that the proposed clustering algorithms stand for the most reliable method for brain image segmentation. For real brain MR image, we classified the image pixels into two classes using the proposed algorithm with Riesz distribution: tumour tissue and healthy tissue. The dark black area corresponds to the tumour with a few pixels are classified wrongly. Obtained experimental results from brain real image indicate excellent agreement and accuracy with Riesz distribution.

This confirmed further by the quantitative segmentation results in terms of the evaluation criteria which are portrayed in Table 5 below.

Table 5. The MCR and PSNR values by varying mixture with brain images.

Figures	Mixture Model	$\mathit{MCR}$	$\mathit{PSNR}$
$(8)$	Mixture Riesz	$0.023$	$29.42$
	Mixture Wishart	$0.042$	$25.56$
$(9)$	Mixture Riesz	$0.002$	$27.54$
	Mixture Wishart	$0.014$	$22.13$

Resting upon Table 5, it is clearly noticed that statistical metrics values of the proposed segmentation approach are much closer to the ideal segmented image. The smallest MCR values as well as the highest PSNR values specify the high quality of segmented image using proposed algorithm with mixture Riesz.

For this reason, we may assert that our approach fulfils a good interpretation and a high accuracy.

Remark 1. As for the Riesz probability distribution, we suggest using a mixture of the inverse Riesz distribution (Tounsi and Zine 2012) for real life applications like image segmentation. Thus, a classification method based on the Expectation-Maximization algorithm within a maximum likelihood framework needs to be developed.

4. Conclusion and discussion

In the current work, we elaborated an iterative algorithm based on the Riesz mixture model and the EM algorithm for settling image segmentation problem. It is an extension to previous research which was confined uniquely to theoretical results of the Riesz distribution without real application. Moreover, our research corroborates the efficiency of the matrix-variate Riesz probability distributions on symmetric matrices in different fields such as image processing. Synthetic and real images are invested to demonstrate the effectiveness of the proposed method, using both visual presentation and numerical evaluation as Misclassification Rate, Peak Signal to Noise Ratio, Adjusted Rand Index, Rand index, Mirkin Metric and Hubert Index. Riesz distribution proved to better characterise and distinguish the regions with different textural information than the Wishart one. Compared to the classical iterative method based on the Wishart distribution, the proposed strategy can effectively improve the computation efficiency as well as segmentation accuracy. With real brain images, we get two classes: class of brain tumour and class of normal tissue. According to experts, the Riesz distribution has a good characterisation which can be invested to characterise the classes in heterogeneous areas. Hence, the proposed method can obtain more accurate and visually pleasing results than the classical statistical methods. The initialisation of the parameters stands for the first step in the proposed segmentation method and represents the most significant stage which is extremely influential in terms of the final estimation results. The unsatisfactory segmentation results can be obtained when the initial parameters are not well established. As a final note, we would assert that in future research, more diligent efforts would be performed in this area so as to overcome this problem.

Acknowledgements

The authors would like to express their gratitude to the reviewers for the time taken to review the manuscript and their valuable comments and suggestions, which have undoubtably improved the quality of the manuscript. The authors express their hearty gratitude to the Radiologist Dr. Moez Bradai (Magnetic Resonance Imaging Center, Road of El-An Km 0.5, 3000 Sfax, Tunisia) for providing useful brain MR images.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional information

Funding

The authors declare that there was no funding received for this research.

Notes on contributors

Mariem Tounsi

Mariem Tounsi is an Assistant Professor in Sfax National Engineering School, University of Sfax, Tunisia. She received Ph.D. degree in Mathematics from Faculty of Sciences of Sfax, University of Sfax, Tunisia, in 2013. She also received M.S. degree in Applied Mathematics from Faculty of Sciences of Sfax, in 2009. Her research interests are related to multivariate analysis, random matrices, statistical inference and theory of probability.

Mouna Zitouni

Mouna Zitouni is an Assistant Professor in Higher Institute of Computer Sciences and Multimedia of Gabes, University of Gabes, Tunisia. She received the Ph.D. degree in Mathematics at Faculty of Sciences of Sfax, University of Sfax, Tunisia, in 2017 and the M.S. degree in Applied Mathematics at Faculty of Sciences of Sfax, in 2013. Her research interests include mixture models and image processing, mainly segmentation and classification techniques.