Improving Parotid Gland Tumor Segmentation and Classification Using Geometric Active Contour Model and Deep Neural Network Framework

1.png

In this work, an advanced classification scheme combined with applicable textural region characteristics is proposed. The first tasks of the proposed method were presented in our previous work [28]. In this contribution, a continuation of the proposed work [13] is presented where widely of information are treated to obtain effective performances in images denoising scheme. The proposed method [28] is included a chief part. Different images pre-processing strategy is first applied by the use of filtering technique to reduce the noise in used images. Second strategy, appreciated parotid gland segmentation is achieved to detect the lesion region based on the GAC. Figure 2 summarizes the overall stages of the proposed analysis scheme.

2.png

Figure 2. The proposed methodology

Hence, the contribution of this work is divided into two principle folds:

•   MR Images segmentation: Appreciated Parotid gland lesion segmentation is achieved to detect the region of interest based on the GAC.

•  Intelligent Classification framework: An image categorization strategy is projected in order to classify the MR subjects into two groups: parotid gland lesion level 1 (PGLL1) and parotid gland lesion level 2 (PGLL2) datasets. The Deep Neural Network (DNN) method is applied for classification step. Here, the previous detection results of PG lesion boundary are employed to characterize the anomaly using T1-MR images [13]. The DNN model classifier is performed to recognize MR images in different cases. This method is primarily recipient for precocious detection of PG defect progression using early period assessment.

2.2.1 Image pre-processing Stage using BM3D technique

In this phase, a filtering method based on BM3D algorithm is implemented to improve the image quality and facilitate subsequent work steps. The main structure of the BM3D algorithm can be split into two major steps. The first step estimates the de-noised image using hard thresholding during the collaborative filtering and can be divided into three sub-steps. These sub-steps start by finding the image patches P(P) similar to given ones and grouping them into a 3D block, then a 3D linear Wavelet transform and shrinkage of coefficient are applied to these blocks. Finally, an inverse 3D transform is applied using the following equation.

$P(P)^{\text {hard }}=\tau_{3 D}^{\text {hard }-1}\left(\gamma\left(\tau_{3 D}^{\text {hard }}(P(P))\right)\right)$    (1)

where, γ is a hard thresholding operator with threshold $\lambda_{3 D}^{\text {hard }} \sigma$ such as σ designates the zero means Gaussian noise variance.

$\gamma(x)= \begin{cases}0 & \text { if }|x| \leq \lambda_{3 D}^{\text {hard }} \sigma \\ x & \text { otherwise }\end{cases}$    (2)

The final 3D block consists of the simultaneously filtered 2D image patches. The last step is the aggregation operation it gives an estimate for each used patch. These estimates are saved in a buffer.

$\forall Q \in P(P), \forall x \in Q,\left\{\begin{array}{l}v(x)=v(x)+w_{P}^{h a r d} \mu_{Q, P}^{h a r d}(x) \\ \delta(x)=\delta(x)+w_{P}^{h a r d}\end{array}\right.$    (3)

where:

- ν and δ designate respectively the numerator and the denominator part of the basic estimate of the image obtained at the end of the grouping step

- $u_{Q, P}^{\text {hard }}$ designate the estimate of the value of the pixel x obtained during collaborative filtering.

$w_{P}^{\text {hard }}= \begin{cases}\left(N_{p}^{\text {hard }}\right)^{-1} & \text { if } N_{p}^{\text {hard }} \geq 1 \\ 1 & \text { otherwise }\end{cases}$    (4)

Such as $N_{P}^{\text {hard }}$ refer to the non-zero coefficients in the 3D block obtained after hard thresholding mentioned above.

The obtained basic estimate after this first step is given by the following equation.

$u^{\text {basic }}(x)=\frac{\sum_{P} w_{P}^{\text {hard }} \times \sum_{Q \in P(P)}  \quad \chi Q(x) u_{Q, P}^{\text {hard }} \quad (x)}{\sum_{P} w_{P}^{h a r d} \quad \times \sum_{Q \in p(P)} \quad \chi Q(x)}$    (5)

This estimate is simply obtained by dividing the two buffers, which are the numerator and the denominator, element by element with respect to the following condition.

$\chi Q(x)= \begin{cases}1 & \text { if and only if } x \in Q \\ 0 & \text { otherwise }\end{cases}$    (6)

The second main step is carried out by using Wiener filtering. This second step imitates the first one with two differences taken into account [2]. The filtered patches are compared instead of the original patches. The new 3D group is processed by Wiener filtering instead of a simple threshold.

Three sub steps are carried out in this stage: designated grouping, collaborative filtering, and aggregation. During the grouping operation, the patch-matching is handled on the basic estimate $u^{\text {basic }}$. In fact, a set of similar patches are collected to form 3D groups.

$P^{\text {basic }}(P)=\left\{Q: d(P, Q) \leq \tau^{\text {wiener }}\right\}$    (7)

where:

- $P^{\text {basic }}(P)$ and P(P) are respectively obtained after the stacking of the similar patches from the basic estimate and the original noisy image.

- $\tau^{\text {wiener }}$ is the distance threshold for d under which two patches are assumed similar.

After the grouping operation, which is achieved by finding the 3D blocks, the collaborative filtering can be started. The filtering operation of P(P) called Wiener collaborative filtering is achieved by using Wiener filtering. In fact, the resulting image is produced by multiplying element by element the 3D transform of the noisy image with the Wiener coefficients defined as follow.

$w_{p}(\xi)=\frac{\left|\tau_{3 D}^{\text {wien }}\left(P^{\text {basic }}(P)\right)(\xi)\right|^{2}}{\left|\tau_{3 D}^{\text {wien }}\left(P^{\text {basic }}(P)\right)(\xi)\right|^{2}+\sigma^{2}}$    (8)

The estimate of the 3D group is obtained by the Eq. (9).

$P^{\text {wien }}(P)=\tau_{3 D}^{\text {wien }^{-1}}\left(w_{P} . \tau_{3 D}^{\text {wien }}(P(P))\right)$    (9)

The last minor step is the aggregation which takes place when the collaborative filtering is achieved. This step aims to store in a buffer the estimate for every pixel.

$\forall Q \in P(P), \forall x \in Q,\left\{\begin{array}{l}v(x)=v(x)+w_{P}^{\text {vien } } u_{Q, P}^{\text {wien }}(x) \\ \delta(x)=\delta(x)+w_{P}^{\text {wien }}\end{array}\right.$    (10)

where:

- τ and δ designate respectively the numerator and the denominator part of the final estimate image which is obtained in the end of the first step.

- $u_{Q, P}^{\text {wien }}(x)$ is the estimation of the value of the pixel x belonging to the patch Q obtained during the collaborative filtering of the reference patch P.

$w_{P}^{\text {wien }}=\left\|w_{P}\right\|_{2}^{-2}$    (11)

The final estimate image is given according to the Eq. (12):

$u^{f_{\text {inal }}}(x)=\frac{\sum_{P} w_{P}^{\text {wien }} \sum_{Q \in P(P)} \quad \chi Q(x) u_{Q, P}^{\text {wien }} \quad (x)}{\sum_{P} w_{P}^{\text {wien }} \sum_{Q \in P(P)} \quad \chi Q(x)}$    (12)

Following the first major step, the final estimate is simply obtained by dividing the two buffers, which are the numerator and the denominator, element by element with respect to the condition mentioned in Eq. (12).

All these steps are summarized in the Figure 3.

3.png

2.2.2 PG lesion boundary detection

The objective of the first stage of the proposed automatic PG lesion segmentation technique is the boundary contour extraction of the lesion region from the remainder (background). This task is based on a geodesic active contour model initialized by an elliptic curve that estimates the lesion region in the PG.

The GAC approach [29] is applied to segment the object edge in a grayscale image by finding a curve C that decrease the following energy function:

$E_{C}(C)=\int_{0}^{L(C)} g(C(q)) d s+v \int_{\Omega} g(C(q)) d A$    (13)

where, $C(q)=(x(q), y(q)$, is a curve in Ω with a position parameter $q \in[0,1]$ and L represents the Euclidean distance. The distance ds is the Euclidean metric, dA is the parameter of area and v is a parameter which controls the area of the active contour. Ωϵ$R^{2}$ is the 2D space of the gray level image, and g is the stopping function of the curve (edge-detector) for gray image. The progress of the curve is controlled by the stopping function represented by the following form:

$g=\frac{1}{1+\left|\nabla\left(G_{\sigma} * I\right)\right|}$    (14)

where, I is the image intensity and G represents the Gaussian kernel. The curve is categorized by two parameters the normal direction N with a convergence parameter F:

where, Gσ is the Gaussian kernel with standard deviation σ. The curve is characterized by two parameters; the normal direction N with speed convergence F:

$\frac{\partial C(t)}{\partial t}=F . \overrightarrow{\mathrm{N}}$    (15)

The level set method, introduced by Osher and Sethian [19] to find a numerical solution of this problem (Eq. (1)), consists on representing the contour C(t) as being the zero level of two dimensional Lipchitz continuous function ϕ(p, t), called the level set function. In other words, the contour C(t) is the set of points P(x, y) such as:

$C(t)=\{p / \phi(p, t)=0\}$    (16)

where, ϕ is the signed distance function that gives the location of a point p with respect to the contour C:

$\phi(p)= \begin{cases}+d, & p \text { is outside } C \\ 0, & p \text { in } C \\ -d & p \text { is inside } C\end{cases}$    (17)

The level set function ϕ for the contour point must have the zero value at each time:

$\phi(C(t), \mathrm{t})=0$    (18)

The derivative of the level set function is given by:

$\frac{\partial \phi(C(t), t)}{\partial t}=\frac{\partial \phi(C(t), t)}{\partial t}+\frac{\partial \phi(C(t), t)}{\partial C} \cdot \frac{\partial C}{\partial t}=0$    (19)

The Eq. (19) is written as:

$\frac{\partial \phi}{\partial t}+\nabla \phi \cdot \frac{\partial C}{\partial t}=0$    (20)

From the Eq. (19), the gradient of ϕ function represents the normal N vector of the contour and the evolution rate is a direction normal to the curve:

$F=\frac{\partial C}{\partial t} . \vec{N}$    (21)

The Eq. (15) is integrated in (20), and, we obtain the evolution equation given by:

$\frac{\partial \phi}{\partial t}+F .\|\nabla \phi\|=0$    (22)

The curve evolution is generally a function of the image and the curve characteristics [30]. The original formulation of the GAC model [19] is given by:

$\frac{\partial \phi}{\partial t}=g(\kappa+v) .\|\nabla \phi\|$    (23)

where, the product g(κ+ν) fixed the evolution speed of level sets of ϕ(p, t) along their normal direction and κ is the curvature which is used here to accomplish the smoothness of the contour.

The curvature is given by [17]:

$\kappa=\operatorname{div}\left(\frac{\nabla \phi}{\|\nabla \phi\|}\right)=\frac{\phi_{x x} \phi_{y}^{2}-2 \phi_{x} \phi_{y} \phi_{x y}+\phi_{y y} \phi_{x}^{2}}{\left(\phi_{x}^{2}+\phi_{y}^{2}\right)^{3 / 2}}$    (24)

The mathematical curve solution of the curve evolution problem is given by the fast level set formulation proposed in [30]:

$\left\{\begin{array}{l}\frac{\partial \phi}{\partial t}=\delta_{\varepsilon}(\phi)\left(\lambda \operatorname{div}\left(g \frac{\nabla \phi}{|\nabla \phi|}\right)+v g\right)+\mu \operatorname{div}((1-1 /|\nabla \phi|) \nabla \phi) \\ \phi(x(t), y(t), 0)=\phi^{0}(x(t), y(t))\end{array}\right.$    (25)

m is a positive parameter that controls the penalization term $\operatorname{div}((1-1 /|\nabla \phi|) \nabla \phi$ which is applied to fix the level set function close to the signed distance function during the curve evolution (see [17]). Where λ presents a constant controlling the curvature term and $\delta_{\varepsilon}(\phi)$ is a parameter that regularizes the Dirac function with a width controlled by e [17]. $\phi^{0}(x, y)$ is the initial contour in each frame of the MR images.

Finally, the discrete level set function is given by:

$\phi^{k+1}(x, y)=\phi^{k}(x, y)+\tau L\left(\phi^{k}(x, y)\right)$    (26)

where, $L\left(\phi^{k}(x, y)\right.$ presents the approximation of the right-hand side in (24) and τ is the time step.

The initial curve is constructed by an ellipse mask selected with five ellipse parameters (x_c, y_c, a, b, θ) where $x_{c}$ and, $y_{c}$ represent the circle centroid, the two ellipse radius (a, b) and θ is to the ellipse angle, which are achieved from the Otsu segmentation method [31]. The idea of this latter method is to detect the PG lesion before the initialization of GAC model. Since the PG presents a low grayscale level, a maximum intensity value is fixed as a threshold based on minimizing the intra-class variance and maximizing the ratio between the interclass variance and the total variance, each frame segmented is composed of two different regions: the region of interest ROI (the Lesion in PG) and the background. Here, to speed up the convergence of the GAC algorithm, we optimized the number of iterations of the curve evolution using the ellipse fitting error. The relative function of error $f_{\text {iter }}$ is computed as follow:

$f_{\text {iter }}=\frac{\text { cover }_{\text {ellip }}-\text { Lesion }_{\text {rough }}}{\text { Lesion }_{\text {rough }}} \times 100$    (27)

where, cover_ellip represents the area that the initial ellipse covers, and the Lesion_rough is the rough lesion area obtained by fixed thresholding. The iteration number noted Iter was computed and reduced as given in the following equations:

iter $=20$ if $\quad f_{\text {iter }} \leq 0.09$    (28)

iter $=30$ if $0.09 \leq f_{\text {iter }} \leq 0.12$    (29)

For the analysis of subsequent MR image frames, the lesion region is chosen based on the current frame segmentation result. In fact, to remove the influence of background artifacts and to limit the computation time of the contour evolution, the initial curve is detected based on the current result.

2.2.3 Classification of PG lesion using DBN-DNN classifier

In this study, the deep neural network method is realized in order to classify the extracted PG lesion boundary and to distinguish between PG lesion levels (level 1 and level 2). Several classification methods can be employed for pattern recognition. Nevertheless, for multiple categorization works [32, 33], it has been confirmed that DNN presents a high-quality tool. To acquire authentic classification results, we have selected the PG region of different images for assembling the training dataset of the network. Consequently, both learning and training neural network stage are performed using the morphological information (texture, area…) of PGs lesion region.

The DNN classification phase is composed into two diverse tasks: (i) the training task that included the greater part of the database (184 subjects randomly chosen using two cases (PGLL1 and PGLL2) which allow the learning and the training of the network. At the same time as (ii) the test task is reserved to validation part that involved the smaller part of the database (the remaining subjects: 122 subjects) taking into account that the examples of the validation set do not belong to the training set. In addition, applying the back-propagation (BP) algorithm, the network is learned with conjugate gradient training method [34, 35].

The neural network structure is composed of two inputs classes, are either (0.95) for (PGLL1), (-0.95) for (PGLL2). In this network, the hyperbolic tangent sigmoid transfer function is applied for all neurons as activation function. The maximum number of iterations (20 epochs of training) is reached or the mean square error (MSE) is decreases to a mean of 8.10^-2.

To improve classification results, all previous works have focused their diagnosis on visual inspection of segmented PG lesion region. The superiority of this study is highlighted by the use of detected ROI by some customized processing tools, and the use of the classification technique using deep learning strategy which was not detained before.

(1) Classification using deep neural network

In this study, the DBN-DNN method is performed to differentiate between two categories of parotid gland lesions. Yet, for numerous classification works [36, 37], it has been proved that DNN represents a high quality methodology for pattern recognition. The deep learning step of the network is performed using PG lesion characteristics. Taking into account the similarity between the characteristics of different PG tumor regions, the DNN procedure is considered using the feature number in the input learning set to obtain the most discriminative features distorted to a suitable illustration of the database used (1×(2 output classes)).

The deep learning step of the network is performed using image feature vectors. Taking into account the resemblance between different characteristics extracted from each PG class in the MR image, the DBN-DNN procedure uses a number of features in the input learning set to obtain the most discriminative characteristics.

The cross-validation experiment is used to give the best DBN-DNN structure. The training-segment vectors sets are divided into five equally-sized sub datasets (five folds). Each iteration contains four folds for the training and one fold for the validation. The tangent sigmoid transfer function is applied for all neurons. Five iterations are achieved for the training and validation phase. After the selection of the hyper parameters, the DBN-DNN classifier was trained and evaluated using the test set (see Figure 4).

4.png

Figure 4. Diagram of the cross-validation and the optimisation process

5.png

Figure 5. The structure of the DBN classifier

(2) Constructing the DBN classifier

The architecture of the DBN classifier is shown in Figure 5, which is composed of three RBMs (RBM1, RBM2 and RBM3) and an output layer which presents the two PG’s lesion levels (PGLL 1 and 2). Firstly, we have to train in an unsupervised way in order to train the DBN classifier. Then, the trained DBN is also combined with the output layer to form the input of the deep neural network (DNN). The DBN is obtained and the DNN is trained by the back-propagation (BP) algorithm in a supervised way.

(3) Training DBN with RBMs

1) Training restricted Boltzmann machine (RBM)

A restricted Boltzmann machine (RBM) is a specific type of random neural network design that has two-layer architecture. The first input layer is called visible layer, where the second layer is called hidden layer. Between two layers, nodes are fully connected, nevertheless in the same layer there is no connection. Though, we constitute a bipartite structure. Figure 6 shows that the bottom layer contains visible variables (nodes) and the top layer contains hidden variables (nodes). The symmetric interaction terms between the visible variables and the hidden variables are represented by the matrix W.

6.png

Figure 6. The structure of the RBM

The energy function of the joint configuration can be defined by:

$E(v, h ; \theta)=-\sum_{i j} W_{i j} \nu_{i j} h_{j}-\sum b_{i} v_{i}-\sum a_{j} h_{j}$    (30)

where, $\theta=\{W, a, b\}$ represents the model parameters, $a_{i}$ is the bias of visible unit i, and $b_{i}$ is the bias of hidden unit j.

The joint probability distribution of a certain configuration is achieved by the Boltzmann distribution (and the energy of this configuration):

$P_{\theta}(v, h)=\frac{1}{Z(\theta)} \exp (-E(v, h ; \theta))$    (31)

where, Z(θ) is the normalization constant.

If ν=(v₁, v₂, …, v_i) is an input vector to the visible layer, the binary state $h_{j}$ of the hidden unit j is fixed to 1 with the probability as follows:

$P\left(h_{j}=1 \mid v\right)=\operatorname{sigmoid}\left(\sum_{i} W_{i j} v_{i}+a_{j}\right)$    (32)

With the states of the hidden units, the binary state $v_{i}$ of visible unit i is set to 1 with the probability below:

$P\left(v_{i}=1 \mid h\right)=\operatorname{sigmoid}\left(\sum_{j} W_{i j} h_{j}+b_{i}\right)$    (33)

The RBM structure is frequently trained as follows:

• Step 1. The states of the visible units are set according to the training data.
•  Step 2. Calculating the binary states of the hidden variables by Eq. (31).
• Step 3. After determining the states of all the hidden units, the states of all visible units are determined
• Step 4. The gradients of W are evaluated by the contrastive divergence (CD) learning algorithm, then the gradient descent algorithm is carrying out to update the parameters W,a,b.