ANN Training Using Granulometric Data from Preprocessed MRI with Morphological Transformations to Detect Multiple Sclerosis Lesions

1. Introduction

Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system (CNS), characterized pathologically by demyelination and clinically by episodes of neurological dysfunction disseminated in space and time, which produces a wide variety of symptoms such as alteration of the sensitivity [1]. Symptoms are determined by the location of demyelinating lesions along the CNS [2]. Initially, MS pathology was defined as an inflammatory process associated with focal plaques of primary demyelination in the white matter of the brain and spinal cord [3]. A reliable and accurate diagnosis of MS is necessary to introduce early treatments for the disease. Disease-modifying drug therapies help to control symptoms and prevent progression [4,5]. Currently, magnetic resonance imaging (MRI) is the most common tool for diagnosing MS, understanding the course of the disease, and analyzing the effects of treatments [6,7].

MRI is a non-invasive method of obtaining detailed images of the internal structure of the body, such as organs and tissues. MRI uses radiofrequency (RF) radiation in the presence of controlled magnetic fields to generate cross-sectional images of the body. MRI is obtained by placing the patient inside a large magnet, which induces a relatively strong external magnetic field. This causes the cores of the body atoms to align with the magnetic field and then the RF signal is applied. The energy released from the body is detected and used to create the MRI by computer [8]. As with any other data acquisition system, in the generation of MRI, there may be a component not correlated with the desired signal, known as noise or random signal. This noise signal is generally caused by spontaneous fluctuations such as the thermal movement of free electrons within real or equivalent electrical components [9]. The signal-to-noise ratio (SNR) is essential for evaluating image quality and determining the use of image processing techniques such as noise removal [10]. On the other hand, MS diagnosis using MRI requires a long time due to the large number of images to be labeled (T1-weighted, T2-weighted, and fluid-attenuated inversion recovery (FLAIR) MRI), so manual diagnostic is susceptible to errors. Also, during the scan the Ghosting effect can happen, which is caused by objects in movement [11]. This can affect image quality and make it difficult to identify areas of CNS injury. Therefore, it is necessary to implement alternative methods of MRI automatic processing to remove noise and other undesirable components.

Most image processing applications require extensive analysis of objects within an image. Segmentation is to divide an image into regions of interest according to their features such as grayscale, color, spatial texture, and geometric shapes [12]. The most common form of segmentation is binary, where each pixel is classified as belonging to the foreground or background. One of the most useful methods for segmenting images is thresholding, which includes choosing an intensity value and then classifying pixels below this value as false and pixels above it as true [13]. On the other hand, mathematical morphology (MM) is a classic image analysis technique that helps to filter the image and extract relevant structures. Granulometry is an MM image measurement tool that determines the size distribution of objects in an image without explicitly segmenting each object [14]. So far, various learning models have been proposed for automatic segmentation of MS lesions in MRI, such as computational algorithms based on convolutional neural networks (CNN) [15]. Ghosh et al. [16] proposed four convolutional encoder networks (CENs) with different network architectures (U-Net, U-Net++, Linknet, and Feature Pyramid Network) to determine the optimal MRI sequence to reduce automatic segmentation times in MS lesions detection. La Rosa et al. [17] proposed a generative adversarial network (GAN) that retrospectively generates realistic magnetization prepared 2 rapid acquisition gradient echoes (MP2RAGE) uniform images (UNI) from MPRAGE images to improve automatic segmentation of MS lesions and tissues. Bandyopadhyay et al. [18] implemented a multilevel thresholding method using evolutionary metaheuristics. The proposal incorporates the concept of altruism into the Harris Hawks Optimization (HHO) algorithm to improve its exploitation capabilities. Macin et al. [19] developed a machine learning (ML) model for the MS diagnosis using brain MRI (axial and sagittal). Features were generated with a fixed-size patch-based (exemplar) feature extraction based on local phase quantization (LPQ). The resulting exemplar multiple parameters LPQ (ExMPLPQ) features were concatenated to form a final feature vector. De Oliveira et al. [20] proposed a deep learning (DL) technique for volumetric quantification of lesions in MRI of MS patients using automatic brain and lesion segmentation by two CNNs. The first CNN was used to perform brain extraction, and the second was for lesion segmentation. Acar et al. [21] proposed a CNN model for identifying MS lesions in brain FLAIR MRI. Hashemi et al. [22] implemented a method to segment MS lesions on FLAIR and T2 MRI using a modified U-Net and modified attention U-Net. Wang et al. [23] proposed a complete CNN U-Net for automatic segmentation of MRI, based on spine features and the contrast between gray levels of intervertebral discs and vertebrae. Rondinella et al. [24] implemented a model that includes a U-Net architecture augmented with a long short-term memory (LSTM) convolutional layer and an attention mechanism, capable of segmenting and quantifying MS lesions detected in MRI. Bose et al. [25] implemented a refined version of the fuzzy c-means (FCM) type 2 technique (EMT2FCM) to isolate diverse tissues in brain MRI, achieved through an improved entropy-based membership function.

While several methods have been proposed for the automatic segmentation of MS lesions in MRI, in some cases only simple data preprocessing, such as image resizing to standardize the input dimensions, has been performed before the algorithm training. Therefore, this paper proposes an alternative MRI preprocessing algorithm that includes two stages:

• Perform morphological opening transformations on brain MRI (MS patients and healthy individuals diagnosed by medical experts), to delete noise and other undesirable components. Also, to compute the granulometry of objects in MRI to characterize the demyelination lesions in the brain white matter caused by MS, and use the resulting data for training two artificial neural networks (ANN) models to predict the MS diagnosis.
• Perform morphological closing transformations on brain MRI (MS patients), to create a reference image (without lesions) and compute the granulometry of objects of an image with lesions and the reference image to be compared. Then, to determine the size of the MS lesions by calculating the differences in granulometry measurements. These measurements could support the decision of specialists to estimate the course or progress of the disease.

2. Materials and Methods

2.1. Database

The analyzed database is the same as [19], and it can be downloaded at: https://www.kaggle.com/datasets/buraktaci/multiple-sclerosis (accessed on 20 October 2023). The acquired dataset includes 200 brain FLAIR MRI (100 axial and 100 sagittal) from MS patients and healthy individuals, who attended the Ozal University Faculty of Medicine in 2021. Medical experts analyzed the sections of axial and sagittal images of MS patients with identifiable MS lesions which were assigned to the MS class, and image sections of healthy individuals with normal appearance, without white matter lesions, and were assigned to the healthy class.

2.2. Mathematical Morphology

MM is a powerful image analysis technique based on set theory, integral geometry, and lattice algebra. The target of morphological operators is to extract relevant structures from an image. It is achieved by testing the image with another set of known shapes called a structuring element (SE). The fundamental morphological operators need the definition of an origin for each SE. This origin allows the positioning of the SE at a given point or pixel: an SE at point x means its origin matches x. The elementary planar isotropic SEs are represented in Figure 1 for octagonal and disk graphs. The pixels that match with the elemental planar SE centered on a given pixel p correspond to the neighbors of p in the graph G plus the pixel itself p, ${\overline{N}}_G\left(p\right)$. The shape and size of the SE are chosen according to the geometric properties of the relevant and irrelevant structures in the image. Irrelevant structures refer to noise or objects to be deleted [26].

2.3. Geodesic Transformations

Some morphological transformations include combinations of an input image with specific SEs. The geodesic transformations approach regards two input images. Initially, a morphological transformation is applied to the first image and then forced to remain either above or below the second image. In this case, morphological transformations are limited to elementary erosions and dilations. In practice, geodetic transformations are iterated until stability [27].

A geodesic dilation involves two images, a marker image and a mask image. By definition, both images must have the same domain and the mask image must be greater than or equal to the marker image. First, the marker image is dilated by the elemental isotropic SE. Then, the resulting dilated image is forced to remain below the mask image. Therefore, the mask image acts as a limit to the spread of the marker image dilation. Let f be the marker image and g be the mask image ($D_f=D_g$ and $f\le g$). The geodesic dilation of size 1 of the marker image f concerning the mask image g is denoted by Equation 1 and it is defined as the minimum point between the mask image and the elemental dilation ${\delta }^{\left(1\right)}$ of the marker image. The geodesic dilation of a binary image is shown in Figure 2.

$${\delta }_g^{\left(1\right)}\left(f\right)={\delta }^{\left(1\right)}\left(f\right)\wedge g.$$

(1)

Geodesic erosion is the double transformation of geodesic dilation concerning the complementation of sets:

$${\epsilon }_g^{\left(1\right)}\left(f\right)=\begin{array}{c}{\left[{\delta }^{\left(1\right)}\left(f^c\right)\wedge g^c\right]}^c\hfill \\ {\left[{\left({\epsilon }^{\left(1\right)}\left(f\right)\right)}^c\wedge g^c\right]}^c\hfill \\ {\epsilon }^{\left(1\right)}\left(f\right)\vee g,\hfill \end{array}$$

(2)

where $f\ge g$ and ${\epsilon }^{\left(1\right)}$ is the elemental erosion. First, the marker image is eroded and then the maximum point is computed with the mask image. Figure 3 shows an example of geodesic erosion. Here, the mask image acts as a limit for the marker image reduction.

2.4. Morphological Reconstruction

In practice, geodesic transformations of a given size are rarely used. However, when they are iterated until stability, very efficient morphological reconstruction algorithms can be created. Geodesic transformations of bounded images always converge after a finite time of iterations (i.e., until the mask image completely obstructs the propagation or reduction of the marker image). The morphological reconstruction of a mask image g from a marker image f is based on this principle. Geodesic transformations by reconstruction allow the deletion of some undesired components without considerably affecting the remaining structures [27].

Like geodesic dilations and erosions, the iteration of the autodual geodesic transformation v always achieves stability for a definition domain of a bounded image. The corresponding transformation is known as the self-dual reconstruction of a mask image g from a marker image f and it is defined as:

$$R_g^v\left(f\right)=v_g^{\left(i\right)}\left(f\right),$$

(3)

where i is such $v_g^{(i+1)}\left(f\right)=v_g^{\left(i\right)}\left(f\right)$. Figure 4 displays self-dual reconstruction. In practice, self-dual reconstruction can be obtained by computing erosion and dilation reconstruction in parallel, the output at a given pixel x alternating between erosion and dilation reconstruction depending on whether $f\left(x\right)$ is above or below $g\left(x\right)$:

$$\left[R_g^v\left(f\right)\right]\left(x\right)=\left\{\begin{array}{cc}\left[R_g^{\delta }\left(f\right)\right]\left(x\right),\hfill & iff\left(x\right)\le g\left(x\right)\hfill \\ \left[R_g^{\epsilon }\left(f\right)\right]\left(x\right),\hfill & otherwise.\hfill \end{array}\right.$$

(4)

2.4.1. Opening and Closing by Reconstruction

For morphological reconstruction, the erosion and dilation of size 1 are iterated until stability is achieved. The geodesic dilation ${\delta }_f^1\left(g\right)$ and the geodesic erosion ${\epsilon }_f^1\left(g\right)$ of size 1 are given by ${\delta }_f^1\left(g\right)=f\wedge \delta \left(g\right)$ with $g\le f$, and ${\epsilon }_f^1\left(g\right)=f\vee \epsilon \left(g\right)$ with $g\ge f$, respectively. When the function g is equal to the morphological erosion or dilation of the original function, the opening by reconstruction ${\tilde{\gamma }}_{\mu B}\left(f\right)\left(x\right)$ is given by:

$${\tilde{\gamma }}_{\mu B}\left(f\right)\left(x\right)=\underset{Untilstability}{\underbrace{{\delta }_f^1{\delta }_f^1\cdots {\delta }_f^1\left({\epsilon }_{\mu B}\left(f\right)\right)\left(x\right),}}$$

(5)

or the closing by reconstruction ${\tilde{\phi }}_{\mu B}\left(f\right)\left(x\right)$ is given by:

$${\tilde{\phi }}_{\mu B}\left(f\right)\left(x\right)=\underset{Untilstability}{\underbrace{{\epsilon }_f^1{\epsilon }_f^1\cdots {\epsilon }_f^1\left({\delta }_{\mu B}\left(f\right)\right)\left(x\right),}}$$

(6)

where B represents the SE and $\mu$ is a size parameter. For this research, a disk-shaped SE with radius $r=\mu$ ($rB=\mu B$) is chosen. For example, a disk of radius $r=5$ represents a disk of 9x9 pixels, therefore, 9 neighbors are analyzed. So, an SE of size n consists of a disk of $(2n-1)\times (2n-1)$ pixels.

2.5. Image Measurements

Image measurements aim to characterize objects in an image by computing some numerical values. For a given criterion, the measurement is discriminant if the obtained values for objects satisfying this criterion differ from those for all other objects. MM provides various image measurement tools such as pattern spectrum or granulometry, direction analysis, texture analysis, shape description, etc. These measures define a vector of features that can be used as input for a statistical analysis or classification method [26].

Granulometry determines the size distribution of objects in an image without explicitly segmenting (detecting) each object, and it is used in some areas to describe the features of the size and shape of individual granules within an object [28]. The formal definition of this concept is presented below.

Granulometry.- Let ${\left({\psi }_{\lambda }\right)}_{\lambda \in Z}$ to be a transformations family depending of a unique positive parameter $\lambda$. This family comprises a granulometry if and only if the following properties are fulfilled: (i) $\forall \lambda$ positive ${\Psi }_{\lambda }$ is increasing, (ii) $\forall \lambda$ positive, ${\Psi }_{\lambda }$ is antiextensional, (iii) $\forall \lambda$ and $\mu$ positives, ${\Psi }_{\lambda }{\Psi }_{\mu }={\Psi }_{max(\lambda ,\mu )}$.

In particular, the family of morphological openings and closings for the numerical instance ${\gamma }_{\mu B}$, ${\phi }_{\mu B}$ with $\mu =\{1,2,...,n\}$ satisfy the previous definition.

The granulometric measurement associated with the clear regions is denoted as $\chi$, and it is computed by:

$$\chi ={\displaystyle \frac{vol\left[{\gamma }_{(\mu -1)B}\left(f\right)\left(x\right)\right]-vol\left[{\gamma }_{\left(\mu \right)B}\left(f\right)\left(x\right)\right]}{vol\left[f\right(x\left)\right]}},$$

(7)

where $vol$ refers to the total number of pixels in an image. The granulometric measurement associated with the dark regions is denoted as $\zeta$ and it is calculated by:

$$\zeta ={\displaystyle \frac{vol\left[{\phi }_{\mu B}\left(f\right)\left(x\right)\right]-vol\left[{\phi }_{(\mu -1)B}\left(f\right)\left(x\right)\right]}{vol\left[f\right(x\left)\right]}}.$$

(8)

2.6. Proposed Algorithm

This paper proposes an alternative MRI preprocessing algorithm including elementary morphological transformations and granulometry computing to characterize MS lesions without affecting the image original dimensions of the image. The proposed algorithm is described below and it is divided into two stages.

• Perform opening morphological transformations on brain images of MS patients and healthy individuals (axial and sagittal MRI), compute the granulometry of objects (Equation 7), and use the resulting data to train two ANN models, applying the following steps:
  • Read the original color image (.png) and convert it to grayscale (e.g. uint8 array 569x1158x3 → uint8 array 569x1158).
  • Perform a morphological opening transformation on the image in grayscale (mask image) to create a marker image using a SE. This operation consists of an erosion followed by a dilation using the same SE. The created SE is disk-shaped with radius r, which matches the geometric properties of the relevant structures of a brain image.
  • Perform an opening by reconstruction transformation on the mask image (Equation 5), using the marker image to identify high-intensity objects in the mask image.
  • Adjust the intensity values of the opened image by reconstruction, which increases the contrast of the output image, to extract relevant structures (MS lesions).
  • Compute the granulometry of objects of the opened image by reconstruction for different radius values ($r=1,2,3,...,15$) of the SE.
  • Enter the granulometry measurements of each sample into two arrays (axial and sagittal) of dimension: 100 samples x 15 features (radius values) = 1500 samples for each array.
  • Train two ANN models with the granulometry measurements, using MATLAB R2023a software, to make predictions of MS diagnosis.
• Perform closing morphological transformations on brain images of MS patients (axial and sagittal MRI) to determine the size of the MS lesions, applying the following steps:
  • Subtract the lesion areas from the mask image (original image) to acquire a brain image without lesions (with holes).
  • Perform a morphological closing transformation on the resulting image (previous step) using a disk-shaped SE with radius r to create a marker image. This operation consists of a dilation followed by erosion using the same SE.
  • Perform a closing by reconstruction transformation on the marker image (Equation 6), using the mask image to fill the holes and create a reference image (without lesions), for making comparisons with the mask image (with lesions).
  • Compute the granulometry of objects of the mask image and the reference image for different values of radius ($r=1,2,3,...,15$) of the SE.
  • Determine the size of MS lesions by computing the differences in granulometry measurements of the mask image and the reference image to support the decision of specialists in estimating the disease progress.

Figure 5 and Figure 6 describe the two-stage procedure of the proposed algorithm. Table 1 and Table 2 present the implemented pseudocode of the preprocessing algorithm.

2.6.1. Artificial Neural Network

The ANN structure consists of an input layer of neurons that receive the sample inputs $X=x_1,x_2,\dots ,x_m$, one or more hidden layers of neurons that convert the values from the previous layer into a weighted linear sum, $w_1{x}_1+w_2{x}_2+\dots +w_m{x}_m$, followed by a not linear activation function used to learn the weights. Then the output layer predicts the class label of the samples [29]. During the learning stage, ANN compares the true class labels with the continuous output values of the nonlinear activation function to compute the classification loss and update the weights.

The dataset is randomly divided into training (80%) and test (20%) sets to validate the ANN models. Then, the test accuracy, the dice similarity coefficient (DSC), true positive rate (TPR) or sensitivity, and true negative rate (TNR) or specificity were computed with the validation dataset to evaluate the performance of prediction models [30,31,32]. DSC is given by,

$$DSC={\displaystyle \frac{2TP}{FP+2TP+FN}},$$

(9)

where TP, TN, FP, and FN represent true positives, true negatives, false positives, and false negatives, respectively. Also, the cross-entropy loss is calculated with different regularization strength values (lambda hyperparameter) to train the ANNs. The weighted cross-entropy loss is,

$$L=\sum _{j=1}^n{\displaystyle \frac{w_{j}log\left(m_j\right)}{Kn}},$$

(10)

where the weights $w_j$ are normalized to sum to n instead 1. The test accuracy is calculated by $Acc=1-L$. The lambda value is adjusted to minimize the loss function. Table 3 describes the configured hyperparameters for the ANNs training.

3. Results

The proposed MRI preprocessing algorithm was implemented on 100 MS sample images (axial) and 100 healthy sample images (sagittal). Also, granulometry was calculated to acquire relevant information about some structures in the brain images. In particular, the morphological opening transformation was used as a filter that satisfies the definition of granulometry and mainly affects the clear regions directly associated with MS lesions.

3.1. Algorithm (Stage 1)

In the first stage, the preprocessing algorithm was implemented in all MS and healthy images (axial and sagittal MRI). Figure 7 and Figure 8 describe the obtained results of the applied opening morphological transformations.

Then, granulometry of objects was computed in all MRI samples for an SE of different radius values ($r=1,2,3,...,15$). Figure 9 compares the granulometry measurements results in brain images of two MS patients (axial and sagittal) and two healthy individuals (axial and sagittal) with similar dimensions respectively (569x1158 pixels).

After granulometry was computed, the resulting data were entered into two arrays (100x15 samples) to train the ANNs. Table 4 and Table 5 present granulometry values of some samples (axial and sagittal). Also, Figure 10 displays the values of the sum of granulometry of 80 MRI samples (axial and sagittal).

Then, the performance metrics (DSC, sensitivity, and specificity) were computed using the confusion matrix results to evaluate the ANN models. Figure 11 shows the counting of predictions of the ANNs.

Also, the cross-entropy loss was calculated with different regularization strength values, the corresponding results are displayed in Figure 12. So, the ANNs were trained using the best lambda regularization strengths. Figure 13 displays the loss function behavior throughout different iterations.

Table 6 compares the performance metrics results of the implemented ANN models.

3.2. Algorithm (Stage 2)

At the second stage of the implemented algorithm, the MS lesion areas (acquired at the intensity adjustment step) were subtracted from the mask image (original image), and the morphological closing transformation of the resulting image (with lesion holes) was computed. Then, the image was closed by reconstruction to fill the holes, and a reference image (without lesions) was created to be compared with the mask image (with lesions). The previous procedure was applied to two MS sample images (axial and sagittal), and the results are displayed in Figure 14 and Figure 15.

Finally, the granulometry of objects of the mask image and the reference image were computed. Figure 16 and Figure 17 show the granulometry results of two MS sample images (axial and sagittal). To determine the size of the MS lesions, the difference of granulometry measurements at each data point was calculated, with a disk-shaped SE for different radius values ($r=1,2,3,...,15$).

4. Discussion

Segmentation is useful for analyzing and identifying damaged tissues or other abnormalities in MRI and for different tasks such as disease diagnosis [12]. Although several computational models have been proposed for the automatic segmentation of MS lesions in MRI, in some cases only simple data preprocessing (image resizing to standardize the input dimensions) has been performed before the algorithm training [14,16,17,18]. Therefore, this paper proposes an MRI preprocessing algorithm based on elementary morphological transformations and granulometry computing. The algorithm helps filter the image by eliminating undesired components without affecting its original dimensions. Also, it improves the characterization of MS lesions by extracting only the relevant structures.

At the first stage of the proposed algorithm, morphological opening transformations were performed on the brain images to identify high-intensity objects (clear objects can be associated with MS lesions), followed by intensity adjustment to increase the contrast and extract the relevant structures. Also, the granulometry of objects in the preprocessed images was computed. The granulometry results in Figure 9 show that the intensity values in MS brains are higher than in healthy brains, so these measurements are directly associated with the lesions. Then, the granulometry values were used to train two ANN models to predict MS diagnosis. The performance metrics results in Table 6 prove granulometry measurements work as efficient predictors to train ANN models. In this way, this algorithm can support the decision of specialists to diagnose MS based on granulometry computing. Also, the loss function behavior in Figure 13 indicates it is not necessary to spend so much time tuning the hyperparameters of a learning model, only by adjusting the regularization strength hyperparameter (lambda) it is possible to minimize the loss. Table 7 compares the performance results of some studies in the analysis of MS lesions in MRI. A limitation of this paper was that some MRIs available in the analyzed database had different dimensions, so only 100 images of MS patients (axial and sagittal) and 100 images of healthy individuals (axial and sagittal) of dimension: 569x1158x3 pixels were chosen.

In the second stage of the implemented algorithm, morphological closing transformations were performed on brain MRI (MS patients), to create an image of a healthy brain approach (reference image). Then, the granulometry of objects of an image with MS lesions and the reference image was computed to determine the size of the lesions. The results in Figure 12 and Figure 13 show it is possible to estimate the size of MS lesions by calculating the difference of granulometry measurements at each data point. These measurements could support the decision of specialists to estimate the progress of the disease and possibly prescribe some specific treatment.

To validate the results of this paper, they were analyzed by a neurologist expert. So, he commented the characterization of the size of the lesions based on granulometry measurements could be useful to monitor a load of activity (new plaques in MS), monitor the progress of other demyelinating diseases such as Devic disease (at the spinal cord level), and monitor tumor growth of gliomas of low grade and even multiform glioblastoma (to evaluate the response to treatment: tumor reduction).