A Hyperspectral Imaging Approach to White Matter Hyperintensities Detection in Brain Magnetic Resonance Images

Abstract

White matter hyperintensities (WMHs) are closely related to various geriatric disorders including cerebrovascular diseases, cardiovascular diseases, dementia, and psychiatric disorders of elderly people, and can be generally detected on T2 weighted (T2W) or fluid attenuation inversion recovery (FLAIR) brain magnetic resonance (MR) images. This paper develops a new approach to detect WMH in MR brain images from a hyperspectral imaging perspective. To take advantage of hyperspectral imaging, a nonlinear band expansion (NBE) process is proposed to expand MR images to a hyperspectral image. It then redesigns the well-known hyperspectral subpixel target detection, called constrained energy minimization (CEM), as an iterative version of CEM (ICEM) for WMH detection. Its idea is to implement CEM iteratively by feeding back Gaussian filtered CEM-detection maps to capture spatial information. To show effectiveness of NBE-ICEM in WMH detection, the lesion segmentation tool (LST), which is an open source toolbox for statistical parametric mapping (SPM), is used for comparative study. For quantitative analysis, the synthetic images in BrainWeb provided by McGill University are used for experiments where our proposed NBE-ICEM performs better than LST in all cases, especially for noisy MR images. As for real images collected by Taichung Veterans General Hospital, the NBE-ICEM also shows its advantages over and superiority to LST.

1. Introduction

White matter hyperintensities (WMHs) are commonly observed on T2W or FLAIR brain MR images of elderly people and related to various geriatric disorders including cerebrovascular diseases, cardiovascular diseases, dementia, and psychiatric disorders [1]. According to [2], WMHs are brain lesions that generally show up as brighter areas and can be visualized by T2W and FLAIR MRI sequences. It is also referred to as Leukoaraiosis and is often found in computed tomography (CT) or MRI of older patients. It is a marker of small-vessel vascular disease. In clinical practice, it is indicative of cognitive and emotional dysfunction, particularly in the ageing population. Its initial discovery was observed in the late 1980s by Hachinski and colleagues [2] who described WMH as patchy low attenuation in the periventricular and deep white matter. Since then, detection of WMH has received considerable interest. Although a supervised method may produce better results, it requires human intervention which is very time-consuming and also suffer the issues of intra- and inter-observer variation [3]. Accordingly, segmentation of WMH has been recently directed to semi-unsupervised and automatic methods which rely on computer assisted tools to help diagnosis to avoid human subjective interpretation. Most importantly, such computer assisted diagnosis can be further used to quantify WMH and calculate its volume [3,4,5,6,7]. However, it also comes with two major issues. One is that most works are based on T1 weighted (T1W), T2W and photon density (PD) or FLAIR images to produce spatial statistics to segment WMH. The other is selection of an appropriate threshold, which ultimately determines the detection results of WMH. Generally, such automatic method is not fully automatic but rather semi-unsupervised because it requires adaptively adjusting threshold values by visual inspection. This paper takes a quite different approach to designing a joint spectral–spatial method that takes advantage of spectral properties provided by MR image sequences to perform subvoxel detection in conjunction with a Gaussian spatial filter to capture spatial contextual information surrounding the WMH detected voxels.

One of the strengths of magnetic resonance imaging (MRI) is its ability in imaging structures of soft tissues. Because an MR image is collected by specifically designed image sequences such as T1W, T2W or PD, it can be considered as a multispectral image [8]. Hyperspectral imaging has recently emerged as an advanced technique in remote sensing to deal with many issues that cannot be resolved by multispectral imaging, specifically, subpixel target detection and mixed pixel classification [9]. Its applications to MRI classification have been also explored in [10,11,12,13,14,15]. However, it seems that using the concept of hyperspectral imaging techniques for WMH detection in brain MRI has not been investigated. This paper presents a new application of hyperspectral imaging in WMH detection of MR brain images.

To expand capability of multispectral imaging to hyperspectral imaging in data analysis, it suffers from insufficient band dimensionality. To resolve this dilemma, a nonlinear band expansion (NBE) process was previously proposed in [16] which used nonlinear functions to produce a new set of nonlinear band images that can be incorporated into the original images to create a new data set. As more such nonlinearly generated images are included, the resulting multispectral image has become a hyperspectral image. In this case, we can take advantage of the well-known hyperspectral subpixel target detection technique, called constrained energy minimization (CEM) [9,17,18,19], to detect the lesion of interest [20]. However, the nonlinearly expanded band images by NBE used in [20] can only capture spectral information nonlinearly but not spatial information. As noted above, effectively detecting the boundary of a lesion may also require spatial information due to the shape of the boundary that is closely related to spatial correlation.

This paper develops a novel NBE approach that expands NBE [20] to produce new band images that can capture not only nonlinear spectral information but also spatial information. Since CEM is a pixel-based technique and does not take into account spatial information. In order for CEM to capture spatial information an iterative version of CEM, to be called Iterative CEM (ICEM), is developed for this purpose. Its idea is to apply a Gaussian filter to a CEM-detection map so that the Gaussian-filtered CEM detection map will contain spatial information to be further fed back as a new band image to create a new image cube. The same process of operating CEM on this new data cube is repeated over again in an iterative manner via feedback loops. To terminate ICEM an automatic stopping rule is also designed, which uses Otsu’s method [21] to threshold the Gaussian-filtered CEM detection map obtained at each iteration as a binary image. If the two consecutive binary images agree within an error threshold measured by Dice similarity index (DSI) [22], then ICEM is terminated and the final Otsu’s thresholded binary image is the desired lesion detection map.

There are several main contributions derived from NBE-ICEM. One is using NBE to create new band images to make a multispectral image a hyperspectral image. Another is including Gaussian filters to capture spatial information. Thirdly, such Gaussian-filtered spatial information is further fed back to be included in the data cube being processed as new images to account for spatial information of detected WMH lesions. Fourthly, the spatial information included in CEM is increased via repeated feedback loop in an iterative manner. That is, the more feedbacks the more spatial information to be included in the data cube for better boundary detection. Fifthly, Otsu’s method is introduced to automatically terminate the iterative process carried out by ICEM. Finally, once the ICEM is terminated, the resulting Otsu’s thresholded binary image is the desired final lesion detection result.

2. Methods

2.1. Nonlinear Band Dimensionality Expansion

An early attempt to expand the original set of band images is to utilize nonlinear functions, for example, auto-correlation and cross-correlation, an idea derived from [16,20]. This type of NBE process is referred to as correlation band expansion process (CBEP). Combining these new CBEP-generated band images with the original set of band images produces a hyperspectral image with sufficient band images.

The CBEP presented in this section is an NBE process using correlation functions to generate new band images from the original set of multispectral images. Its original idea was developed in [16,20].

Correlation Band Expansion Process (CBEP)

Step 1. First-order band image: ${\left\{B_l\right\}}_{l=1}^L$ = set of original band images

Step 2. Second-order correlated band images:

(i)

${\left\{B_l^2\right\}}_{l=1}^L$ = set of auto-correlated band images

(ii)

${\left\{B_k{B}_l\right\}}_{k=1,l=1,k\ne l}^{L,L}$ = set of cross-correlated band images

Step 3. Third order correlated band images

(i)

${\left\{B_l^3\right\}}_{l=1}^L$ = set of auto-correlated band images

(ii)

${\left\{B_k^2{B}_l\right\}}_{k=1,l=1,l\ne k}^{L,L}$ = set of two cross-correlated band images

(iii)

${\left\{B_k{B}_l{B}_m\right\}}_{k=1,l=1,m=1,k\ne l\ne m}^{L,L,L}$ = set of three cross-correlated band images

Step 4. Other nonlinear correlated band images

(i)

${\left\{\sqrt{B_l}\right\}}_{l=1}^L$= set of band images stretched out by the square-root.

(ii)

${\left\{\log \left(B_l\right)\right\}}_{l=1}^L$ = set of band images stretched out by the logarithmic function.

It should be noted that, according to the nonlinear functions described in Steps (1)–(4), the band images generated by CBEP contain only nonlinear spectral information but not spatial information. In what follows, we develop an iterative CEM (ICEM) to address this issue where spatial information can be captured by using a Gaussian filter and feed it back to expand images currently being processed to create a new set of image data cubes.

2.2. Iterative CEM

ICEM, presented in this section, is implemented in conjunction with CBEP in an iterative manner. More specifically, it utilizes CBEP to create new band images via an NBE process. Once CBEP process is completed, a new set of image data cubes is generated for CEM to perform subpixel target detection. To obtain class spatial information, a Gaussian filter is introduced in the CEM-detected maps so that spatial contextual information of data sample vectors can be captured by a Gaussian filter. The resulting Gaussian-filtered CEM-detection abundance fractional map is fed back to create a new band incorporated into NBE to form a new hyperspectral cube which will be further used for re-processing CEM again. The same process is repeated over and over again until a stopping rule is satisfied. This repeated implementation of CEM via feedback loops in an iterative fashion is ICEM.

Specifically, at each iteration, say kth iteration, a Gaussian filter is used to blur $|B{|}_{\mathrm{CEM}}^{\left(k\right)}$ which is the absolute value of CEM-detection abundance fractional map, $B_{\mathrm{CEM}}^{(k)}$. This Gaussian-filtered band image, $|B{|}_{\mathrm{GF}\left(\mathrm{CEM}\right)}^{\left(k\right)}$ provides spatial classification information as similar filters used in [23] and will be further fed back to ${\mathsf{\Omega }}_{\mathrm{NBE}}^{(k)}$ to create a new set of hyperspectral images, ${\mathsf{\Omega }}_{\mathrm{NBE}}^{(k+1)}={\mathsf{\Omega }}_{\mathrm{NBE}}^{(k)}\cup \left\{|B{|}_{\mathrm{GF}\left(\mathrm{CEM}\right)}^{\left(k\right)}\right\}$ to be used by CEM again for next iteration. The same procedure is continued. To terminate the process, an automatic stopping rule is designed. It applies Otsu’s method [21] to $|B{|}_{\mathrm{GF}\left(\mathrm{CEM}\right)}^{\left(k\right)}$ to produce a binary classification map, $B_{\mathrm{binary}}^{(k)}$ that will be used to calculate DSI [22]. If two consecutive DSI values are within an error threshold, ICEM will be terminated and $|B{|}_{\mathrm{CEM}}^{\left(k\right)}$ and $B_{\mathrm{binary}}^{(k)}$ will be the desired final real-valued WMH lesion detection map for visual inspection and binary value detection maps of WMH lesions for quantitative analysis.

In the following, we describe detailed step-by-step implementation of ICEM in great detail.

ICEM

• Initial condition: Let ${\left\{B_l\right\}}_{l=1}^L$ be the original set of band images.
• Use an NBE process to create a new set of nonlinear band images, ${\left\{B_i^{\mathrm{NB}}\right\}}_{i=1}^{n_{\mathrm{NB}}}$ where n_NB is the number of new band images by the NBE process.
• Form a new set of band images, ${\mathsf{\Omega }}^{(0)}={\left\{B_l\right\}}_{l=1}^L\cup {\left\{B_i^{\mathrm{NB}}\right\}}_{i=1}^{n_{\mathrm{NB}}}$. Let $d^{(0)}={\left(d_1,\cdots ,d_L,d_1^{\mathrm{NB}},\cdots ,d_{n_{\mathrm{NB}}}^{\mathrm{NB}}\right)}^T$ be the desired target pixels in ${\mathsf{\Omega }}^{(0)}$. Let ${\delta }_0^{\mathrm{CEM}}$ be CEM using d⁽⁰⁾ and R⁽⁰⁾ which are obtained from ${\mathsf{\Omega }}^{(0)}$. Let $k=1$.
• At the kth iteration, update $d^{(k)}$ and $R^{(k)}={{\displaystyle {\sum }_{i=1}^N{r}_i^{(k)}\left(r_i^{(k)}\right)}}^T$ from ${\mathsf{\Omega }}^{(k)}$.
• Use new generated d^(k) and R^(k) for ${\delta }_k^{\mathrm{CEM}}$ to be implemented on ${\mathsf{\Omega }}^{(k)}$. Let $B_{\mathrm{CEM}}^{(k)}$ be the detection abundance fractional map produced by ${\delta }_k^{\mathrm{CEM}}$.
• Use a Gaussian filter to blur $|B{|}_{\mathrm{CEM}}^{(k)}$ where $|B{|}_{\mathrm{CEM}}^{(k)}$ is the absolute value of $B_{\mathrm{CEM}}^{(k)}$. The resulting image is denoted by Gaussian-filter $|B{|}_{\mathrm{GFCEM}}^{(k)}$.
• Check if $|B{|}_k^{\mathrm{CEM}}$ satisfies a given stopping rule. If no, continue. Otherwise, go to Step 9.
• Form ${\mathsf{\Omega }}^{(k+1)}={\mathsf{\Omega }}^{(k)}\cup \left\{|B{|}_{\mathrm{GFCEM}}^{(k)}\right\}$. Let $k\leftarrow k+1$ and go to Step 4.
• $B_{\mathrm{CEM}}^{(k)}$ is the desired detection abundance fractional map and ICEM is terminated.

Figure 1 delineates how ICEM is processed as a detector where ICEM uses Gaussian filters to smooth CEM-detection abundance fractional maps and feeds back Gaussian-filtered CEM detection abundance fractional maps to provide spatial information for re-processing CEM iteratively. By gradually increasing more spatial information through feedback loops the boundaries of WMH lesions can be detected more effectively. It should also be noted that, if a particular NBE technique is used such as CBEP, then NBE-ICEM can be specified by CBEP-ICEM.

2.3. Stopping Rule for ICEM

To effectively terminate ICEM, DSI defined in [22] as

$${\mathrm{DSI}}^{(k)}=\frac{2|S_k\cap S_{k-1}|}{|S_k\cup S_{k-1}|}$$

(1)

is used a stopping criterion where |S| is size of a set S, S_k and S_k−1 are the kth thresholded binary image of the kth CEM detection abundance fractional map, $|B{|}_k^{\mathrm{CEM}}$ and k − 1st thresholded binary image of the k − 1st CEM detection abundance fractional map, $|B{|}_{k-1}^{\mathrm{CEM}}$. Figure 2 depicts a flow chart of a stopping rule using DSI with $\epsilon$ as a prescribed error threshold.

2.4. Algorithm for NBE-ICEM

Using Figure 1 and Figure 2, an algorithm developed to implement ICEM in conjunction with NBE can be described as follows. Figure 3 describes a graphic flow chart of implementing NBE-ICEM.

NBE-ICEM

• Initial conditions:
  • For each class, find its sample mean to calculate the desired signature d for the particular class.
  • Select the values of the parameter σ used for Gaussian filters in ICEM,
  • Prescribe an error threshold $\epsilon$ for DSI in Equation (1)
• Use the NBE process described in Section 2.1 to generate a set of nonlinear band images, ${\left\{B_l^{\mathrm{NB}}\right\}}_{l=1}^{n_{\mathrm{NB}}}$.
• Apply ICEM described in Figure 1 to ${\mathsf{\Omega }}^{(0)}={\left\{B_l\right\}}_{l=1}^L\cup {\left\{B_l^{\mathrm{NB}}\right\}}_{l=1}^{n_{\mathrm{NB}}}$.
• Use DSI described in Figure 2 as a stopping rule to terminate ICEM.
• Output $|B{|}_{\mathrm{CEM}}^{\left(k\right)}$, which is real-valued, and $B_{\mathrm{binary}}^{(k)}$, which is binary-valued, to produce a confusion matrix for classification.

3. Results

3.1. Synthetic Image Experiments

To conduct an objective quantitative study, the synthetic MR brain images containing multiple sclerosis (MS) lesions obtained from the MR imaging simulator of McGill University, Montreal, Canada were used for experiments [24]. MS lesions are typically hyperintense on T2W or FLAIR sequence image. Figure 4a–c shows a slice MR brain image along with the ground truth of MS lesion shown in Figure 4d. The MR brain images are acquired by the modalities of T1W, T2W and PD with specifications provided in BrainWeb site [24]. The thickness of slice is 1 mm with size of $181\times 217\times 181$. Each slice is specified by INU (intensity non-uniformity) 0% or 20%, denoted by rf0 and rf20 with six different levels of noise, 0%, 1%, 3%, 5%, 7% and 9%. The noise in the background of the simulated images is simulated by Rayleigh statistics and signal regions are simulated by Rician statistics. The “percentage (%) of noise” represents the ratio of the standard deviation of the white Gaussian noise to the signal for a reference tissue [24] in terms of %. There were 23 MR images from 91 to 113 slices for our study. To implement ICEM, we need to know the desired target signature d. Two ways were suggested to select training samples to calculate d. One is called all slices-selected training samples, which selects a small set of training samples from all MR image slices. The other is called single slice-selected training samples, which selects a small set of training samples from a particular single MR image slice that can be further used to find training samples for entire MR image slices. Its idea was derived from the extrapolation process used in volume sphering analysis (VSA) developed in [14,15].

Table 1. Specifications of parameters used by NBE-ICEM for BrainWeb images.

Band images	T1W, T2W, PD (3 bands)
Correlation Band Expansion Process (CBEP)	3rd order correlated band images
d	found by all slices-selected or single slice-selected training samples
Gaussian window size	$3\times 3$($5\times 5$)
σ used in Gaussian filter	0.1 with window size $3\times 3$ (0.5 with window size $5\times 5$)
Thresholding method	Otsu’s method
error threshold (DSI)	0.80

specifies the values of parameters used for experiments where two Gaussian filters using window sizes of $3\times 3$ and $5\times 5$, and two different σ = 0.1 and 0.5. The experiments were conducted for all MR image slices according to Table 1 where the results obtained by Gaussian window of $5\times 5$ with σ = 0.5 are tabulated in parentheses.

To further evaluate the performance of the proposed NBE-ICEM, a commonly used segmentation approach, called lesion segmentation tool (LST) [25,26], was used for comparative study. It was originally developed for the segmentation of MS lesions and has also been proven to be useful for the segmentation of brain lesions.

Table 2. Averaged DSI values of lesions detection by CBEP-ICEM1, CBEP-ICEM2 and LST over MR image slices 91–113.

	Methods	CBEP-ICEM1	CBEP-ICEM2	LST
Noise/INU Level
n0/rf0		0.865	0.808	0.739
n1/rf0		0.886	0.864	0.749
n3/rf0		0.893	0.863	0.750
n5/rf0		0.806	0.839	0.731
n7/rf0		0.652	0.822	0.693
n9/rf0		0.579	0.801	0.714
n0/rf20		0.861	0.829	0.733
n1/rf20		0.867	0.834	0.753
n3/rf20		0.881	0.827	0.746
n5/rf20		0.814	0.831	0.732
n7/rf20		0.714	0.825	0.694
n9/rf20		0.540	0.806	0.655

tabulates DSI values calculated by Equation (1) averaged over 23 MR image slices 91–113 of lesion detection produced by CBEP-ICEM1, CBEP-ICEM2 and LST for six different noise levels and two INU levels where all slices-selected training samples were used to find the knowledge of d. The shaded DSI values in Table 2 are the best results. As we can see from the table, CBEP-ICEM1 performed better than CBEP-ICEM2 when noise level is low. However, when noise level is high, CBEP-ICEM2 performed better than CBEP-ICEM1. Nonetheless, both NBE-ICEM-based methods, i.e., CBEP-ICEM1 and CBEP-ICEM2, performed better than LST. It should be also noted that, since LST produced real-valued gray scale images, it required a threshold value to segment WMH lesions. The LST results in Table 2 were obtained by manually adjusting threshold values in order to yield the highest detection rate.

Similarly,

Table 3. Averaged DSI values of lesions detection by CBEP-ICEM1, CBEP-ICEM2 and LST over MR image slices 91–113 using slice 102 to select training samples.

	Methods	CBEP-ICEM1	CBEP-ICEM2	LST
Noise/INU Level
n0/rf0		0.798	0.784	0.739
n1/rf0		0.848	0.847	0.749
n3/rf0		0.871	0.858	0.750
n5/rf0		0.776	0.836	0.731
n7/rf0		0.625	0.816	0.693
n9/rf0		0.389	0.778	0.714
n0/rf20		0.844	0.834	0.733
n1/rf20		0.859	0.837	0.753
n3/rf20		0.854	0.814	0.746
n5/rf20		0.811	0.819	0.732
n7/rf20		0.710	0.804	0.694
n9/rf20		0.549	0.799	0.655

also tabulates DSI values calculated by Equation (1) averaged over 23 MR image slices 91–113 of lesion detection produced by CBEP-ICEM1, CBEP-ICEM2 and LST for six different noise levels and two INU levels where single slice-selected training samples were used to find the knowledge of d and the slice 102 was chosen as the desired single slice. The selection of slice 102 is empirical as long as it includes sufficient tissue information, in which case the middle MR image slice can serve as this purpose. The same conclusions drawn from Table 2 were also valid for Table 3, even though the results in Table 3 were slightly degraded compared to the results in Table 2 because the single slice-selected training samples were used. It should be noted that the results of LST in Table 2 and Table 3 were the same because LST did not allow users to select training samples. This disadvantage is further offset by a need of finding an appropriate threshold value to segment lesion out from the background.

For an illustrative purpose, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 only show detection results of WMH lesions of the 97th MR image slice with six levels of noise and 0% INU by two versions of CBEP-ICEM, using the Gaussian window size of $3\times 3$ specified by σ = 0.1 and the Gaussian window size of $5\times 5$ specified by σ = 0.5, referred to as CBEP-ICEM1 and CBEP-ICEM2, respectively, where two sets of training samples selected by all slices and the single 102nd slice were used to calculate the desired target signatures d to implement NBE-ICEM. As we can see by visual inspection against the ground truth in Figure 4d, CBEP-ICEM1 and CBEP-ICEM2 using two sets of training samples, i.e., all slices-selected and single slice-selected training samples, produced very close results and they both also performed better lesion detection than LST did.

Two comments are noteworthy.

• Despite the fact that the training samples used for our proposed NBE-ICEM were selected based on 2D images such as by all slices and a single slice, these training samples were either stacked as voxels from all slices or extrapolated from a single slice as voxels by VSA. Accordingly, NBE-ICEM is actually run on 3D images as image cubes.
• There is an issue in implementing LST. Since it is packaged as a software algorithm, there is no flexibility for users to choose parameters at their discretion. Besides, it cannot implement T1W, T2W or FLAIR images alone. Instead, it must require T1W images as reference images to segment WMHs [26]. Most importantly, it produces real valued gray level images, which require users selecting a threshold value from a range from 0.05 to 0.95 with a step size of 0.05 to detect WMHs. In [26], this threshold value was suggested between 0.25 and 0.4. However, in practical applications, the best value is generally selected manually. Thus, technically speaking, LST is not fully automatic. Specifically, when synthetic images from the BainWeb were used for experiments, it was found that using both T1W and T2W could not segment WMHs. It must use T1W and PD to detect WMHs and the threshold value must be set to around 0.2 to segment WMHs.

3.2. Real Image Experiments

Real MRI brain images were acquired at the Taichung Veterans General Hospital (TCVGH) by Siemens Magnetom Aera 1.5 Tesla (Erlangen, Germany) MR scanner with a 16-channel phase-array head coil. MR imaging protocol included T1W with 3D MPRAGE, T2W and FLAIR. Since T1W, T2W and FLAIR images used for experiments were collected by 3D high resolution sequences with each voxel of size, 1 × 1 × 1 mm³, the interpolation artifacts and partial volume do not have much effect on imaging. However, as a part of trade-off, this also requires additional 2 min for image acquisition. Other imaging parameters used for data acquisition were voxel size of 1 × 1 × 1 mm³, matrix size = 256 × 256 × 176, NEX = 1. According to a clinical visual inspection criterion [27], the WMH lesions can be graded by Fazekas with three grades of Fazekas shown in Figure 11 for illustration.

A total of 111 cases were collected and all the participants have been well-informed and signed their consents. In addition, the study conducted in this paper was approved by the Ethics Committee of Clinical Research, Taichung Veterans General Hospital (IRB number: CE16138A). Among all the 111 cases there are 58 cases of Fazekas grade 1, 44 cases of Fazekas grade 2 and 9 cases of Fazekas grade 3. Thus, in this study, we selected 10 cases from Fazekas grade 1, 11 cases from Fazekas grade 2, and 9 cases from Fazekas grade 3.

As demonstrated by synthetic image experiments, CBEP-ICEM2 was shown to be a better WMH detection technique. Thus, CBEP-ICEM 2 was used in real image experiments.

Table 4. Parameters used by CBEP-ICEM2.

Band	T1W, T2W, FLAIR (3 bands)
CBEP	3rd order correlated band images
d	found by all slices-selected or single slice-selected training samples
Fazekas grade	1	2	3
Gaussian window size	5 × 5
σ used in Gaussian filter	0.5 with window size 5 × 5
Thresholding method	Otsu’s method
stopping threshold (DSI)	0.80

tabulates the values of parameters used by CBEP-ICEM2 where two sets of training samples selected by all slices and the single 90th slice were selected to calculate the desired target signature d to implement NBE-ICEM.

Figure 12, Figure 13 and Figure 14 show the WMH lesion detection results produced by CBEP-ICEM2 and LST for three Fazekas grades, respectively, where Figure 12a, Figure 13a and Figure 14a are original T1W, T2W and FLAIR MR images; Figure 12b, Figure 13b and Figure 14b are iterative WMH lesion detection images by CBEP-ICEM2 along with final WMH lesion detection by Otsu’s method; and Figure 12c, Figure 13c and Figure 14c show comparisons between lesion detections by CBEP-ICEM2 and LST.

As demonstrated in Figure 12, Figure 13 and Figure 14, our proposed NBE-ICEM using two sets of training samples performed very similarly and also better than LST according to clinical visual evaluation criterion, Fazekas grades [27].

4. Discussion

This paper is believed to be the first work ever reported in the literature to attempt to use a hyperspectral subpixel detection, NBE-ICEM, to detect WMHs on MRI. As demonstrated in Table 2 and Table 3 and Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 and Figure 12, Figure 13 and Figure 14 the synthetic and real image experiments confirm significant improvements using NBE-ICEM over the LST method in WMH detection.

In comparison between CBEP-CEM1 and CBEP-CEM2, we found from Table 2 and Table 3 in the synthetic image experiments that, when the noise level is low (0%, 1%, 3%), CBEP-ICEM1 using a smaller Gaussian window performed better than CBEP-ICEM2 using a larger Gaussian window. However, when the noise level is high (5%, 7%, 9%), the conclusion is reversed, i.e., CBEP-ICEM2 performed better than CBEP-ICEM1. It is also interesting to note that CBEP-ICEM1 performed very poorly when noise level reached 7% and above and even worse than LST. Table 2 and Table 3 also shown that it was noise not INU that had an impact on lesion detection. On the other hand, CBEP-ICEM2 generally performed well regardless of noise level if DSI value was set to at least or above 0.8 compared to LST whose DSI values did not go beyond 0.8. The synthetic image experiments suggested that CBEP-ICEM2 was a better technique due to its robustness to noise level and ability in WMH detection.

In addition, based on the results of real image experiments from Figure 12, Figure 13 and Figure 14, there are three interesting findings. Firstly, the number of iterations carried out by CBEP-ICEM is always two for all three Fazekas grades. Secondly, in Figure 12c, CBEP-ICEM2 and LST performed similarly but quite different in Figure 13c and Figure 14c, where the areas of lesions detected by LST were much smaller than CBEP-ICEM2. Thirdly, the iterative images produced in Figure 12b, Figure 13b and Figure 14b by CBEP-ICEM2 showed that including spatial information captured by Gaussian-filtered CEM detection images did improve lesion detection, particularly edge and boundary pixels.

This paper makes several main contributions to WMH lesions detection in MR brain images. First, it develops NBE to resolve two issues arising in WMH detection, insufficient spectral dimensionality and linear non-separability problem. NBE plays a similar role that kernels play in pattern classification such as support vector machine (SVM). Second, it introduces Gaussian filters to be included in CEM to expand capability of CEM in capturing spatial information surrounding CEM-detected WMH lesions. Third, the real-valued CEM-detection abundance fractional maps provide soft decisions for visual inspection. Fourth, Otsu’s method is incorporated in ICEM to produce thresholded binary maps as hard decisions that show WMH lesions detection. This resolves the main issue encountered in LST. Fifth, the feedbacks of Gaussian filtered CEM detection abundance fractional maps allow CEM to perform better detection in WMH lesions when spatial information of lesions is crucial, specifically, their boundaries. Finally, an automatic stopping rule is particularly designed to determine how much spatial information is needed for CEM to perform its best in detection of WMH lesions.

5. Conclusions

In conclusion, this paper develops a novel approach, called NBE-ICEM, for WMH lesions detection in MR brain images. It is derived from a hyperspectral imaging-based subpixel target detection method (CEM), but is rather different from CEM in two aspects. One is an introduction of NBE into CEM, which expands the original MR images by including nonlinearly correlated band images generated by NBE to make a multispectral MR image into a hyperspectral MR image since CEM is a hyperspectral imaging technique. The other is the development of an iterative version of CEM, ICEM, which can feed back spatial information captured by Gaussian filters in an iterative manner. More specifically, it applies a Gaussian filter to CEM-detection maps to produce Gaussian-filtered CEM-detected abundance fractional maps that can be further fed back iteratively to form a new set of MR image cubes which will be used a new data set to be re-processed by CEM to improve WMH lesions detection performance. ICEM can be considered as a joint spectral–spatial filter. The more iterations carried out by ICEM, the more spatial information captured. The work on NBE-ICEM presents a potential and promising technique for WMH lesions detection. It is our belief that there would be new applications of NBE-ICEM to MRI yet to explore in the future, specifically, partial volume estimation for specific tissues of interest.