Skip to main content

Semi-automated brain tumor segmentation on multi-parametric MRI using regularized non-negative matrix factorization



Segmentation of gliomas in multi-parametric (MP-)MR images is challenging due to their heterogeneous nature in terms of size, appearance and location. Manual tumor segmentation is a time-consuming task and clinical practice would benefit from (semi-) automated segmentation of the different tumor compartments.


We present a semi-automated framework for brain tumor segmentation based on non-negative matrix factorization (NMF) that does not require prior training of the method. L1-regularization is incorporated into the NMF objective function to promote spatial consistency and sparseness of the tissue abundance maps. The pathological sources are initialized through user-defined voxel selection. Knowledge about the spatial location of the selected voxels is combined with tissue adjacency constraints in a post-processing step to enhance segmentation quality. The method is applied to an MP-MRI dataset of 21 high-grade glioma patients, including conventional, perfusion-weighted and diffusion-weighted MRI. To assess the effect of using MP-MRI data and the L1-regularization term, analyses are also run using only conventional MRI and without L1-regularization. Robustness against user input variability is verified by considering the statistical distribution of the segmentation results when repeatedly analyzing each patient’s dataset with a different set of random seeding points.


Using L1-regularized semi-automated NMF segmentation, mean Dice-scores of 65%, 74 and 80% are found for active tumor, the tumor core and the whole tumor region. Mean Hausdorff distances of 6.1 mm, 7.4 mm and 8.2 mm are found for active tumor, the tumor core and the whole tumor region. Lower Dice-scores and higher Hausdorff distances are found without L1-regularization and when only considering conventional MRI data.


Based on the mean Dice-scores and Hausdorff distances, segmentation results are competitive with state-of-the-art in literature. Robust results were found for most patients, although careful voxel selection is mandatory to avoid sub-optimal segmentation.

Peer Review reports


High-grade gliomas (HGGs) account for 80% of all malignant primary brain tumors [1]. Standard treatment of HGG consists of complete or partial resection of the gross tumor volume, followed by radiotherapy and/or chemotherapy to attack the remaining tumor cells and to counteract tumor recurrence. So far, recent advancements in glioma research have only had little effect on patient outcome. The 5 year survival rate of the most common histopathological subtypes, anaplastic astrocytoma and glioblastoma multiforme (GBM) are 26 and 5%, respectively [1].

MRI has become the imaging modality of choice for evaluating tumor progression and the efficacy of a chosen treatment strategy. GBM, the most malignant and most common type of glioma, tends to invade the healthy tissue, such that tumor margins extend beyond the imageable component of the tumor based on conventional MRI (cMRI). Numerous recent studies have recommended incorporating additional imaging biomarkers from advanced MRI modalities such as perfusion-weighted imaging (PWI) and diffusion-weighted imaging (DWI) [24]. PWI is used for studying tumor angiogenesis, as HGGs are known to stimulate vascular ingrowth [5]. PWI measurements in the peritumoral region have been reported to differentiate GBM from metastases, suggesting the detection of tumor infiltration [6]. PWI was also found useful in differentiating tumor recurrence from radiation necrosis [7]. DWI assesses the mobility of water molecules within tissue due to Brownian motion. Because glioma infiltration disrupts the organization of the white matter tracts, DWI is potentially useful to characterize the extent of infiltration. Diffusion tensor imaging has been reported to better delineate tumor margins in gliomas than cMRI [8].

Assessment of tumor extent plays a key role at all stages of the treatment process. An outline of the gross tumor volume is made upon planning surgical resection and for defining the radiotherapy target volume. During follow-up, the response assessment in neuro-oncology (RANO) criteria are applied for monitoring tumor growth, assessing the maximal diameter of the lesion in 2 orthogonal directions on axial slices [9]. Clinical practice would benefit from volumetric measurements of the tumor and its subcompartments. With the recent emergence of the dose-painting concept in radiotherapy, segmentation of the active tumor region would allow focusing the radiation energy to the most active part of the tumor in a non-uniform way [10]. Volume measures of the active tumor and necrotic region have been reported to be significant predictors of patient outcome to treatment [11, 12]. Furthermore, the RANO criteria acknowledge that volumetric measurements might be favorable compared to cross product diameters in cases of irregularly shaped tumors, multi-focal tumors and tumors with cystic or necrotic components. However, segmentation of HGGs in MRI is a challenging task, due to their heterogeneous nature: several stages of the disease can occur throughout the same lesion and diffuse boundaries exist between active tumor, necrosis, edema and the surrounding healthy brain. Manual segmentation of 3D images is also time-consuming, which is the main reason why volumetric tumor delineation has not become widespread in clinical practice.

In recent years, significant advancements have been made in the field of automated brain tumor segmentation. Both semi-automatic and fully automatic methods have been proposed [13]. A popular approach is to combine imaging biomarkers from different MRI sequences on a voxel-wise basis, thereby increasing specificity and reducing overlap between tissue classes. Nowadays, supervised classification methods are receiving most attention [14]. These methods rely on an extensive set of training images with manually annotated ground truth to learn decision boundaries between the tissue classes in feature space. The most popular methods include random forests [15, 16], support vector machines [17, 18] and neural networks [19]. Additional constraints may be imposed to further enhance performance, such as spatial consistency of the tissue regions [20, 21].

Unsupervised classification methods are also being considered for tumor segmentation. These methods are very flexible, as they don’t require an extensive training dataset with a uniform acquisition protocol. They learn classification rules directly from the imaging data at hand, based on some similarity criterion. Popular approaches include fuzzy C-means clustering (FCM) [22, 23], Gaussian mixture modeling [24, 25], hidden Markov Random Fields [26, 27] and non-negative matrix factorization (NMF) [28, 29]. Due to the absence of training data, unsupervised methods rely more strongly on the incorporation of prior knowledge to obtain competitive results. As gliomas exhibit a wide variability in terms of size, shape, location and appearance, imposing proper prior knowledge is difficult. Several studies make use of a normal brain atlas to differentiate the pathological region from the healthy tissue structures [25, 30]. This approach assumes that the healthy brain structures are not altered by the tumor, which might not be valid for large tumors deforming the surrounding healthy tissue. Some studies assume the active tumor region to be either enhancing [30] or non-enhancing [22], or hyper-perfused [23]. These assumptions do not hold in general in heterogeneous gliomas with varying degrees of contrast enhancement and/or perfusion. Another approach is to detect tumors by looking at the dissimilarity across the hemispheres, thereby supposing that the tumor volume is restricted to one of the hemispheres [31]. For the removal of false positives, several segmentation algorithms only retain the largest connected component in the segmentation mask, thereby assuming only one tumor volume [22, 32]. Evidently, many unsupervised segmentation studies incorporate prior knowledge in the form of simplifying rules, limiting their applicability to a subset of gliomas with specific characteristics. Some methods are limited to a particular set of MRI parameters, as they rely on one or more specific imaging biomarkers.

Semi-automated segmentation algorithms incorporate prior knowledge in the form of user-specific input, either as an initialization or as a post-processing step. Competitive performance has been reported for semi-automated algorithms [14]. Kwon et al. combined a normal brain atlas with a tumor growth model [33]. User-defined seeding points are used to initiate the tumor growth model, allowing to model multi-focal gliomas as well. Hamamci et al. initialized a cellular automata algorithm based on the maximal tumor diameter as drawn by the user on T1C images [34]. An active level-set surface is then initialized from the user-defined maximal diameter, to impose spatial smoothness on the contour of the pathological region. Havaei et al. initialized a k-nearest neighbors (kNN) classifier with user-selected voxels in the tumor subcompartments as well as in the healthy brain regions [35]. Spatial coordinates of the seeding points were also exploited in the feature set. The main drawback of semi-automated methods is that they don’t provide reproducible results, as the segmentation depends on subjective user input. Robustness against user input variability is therefore a key aspect of these methods. In the current study, we propose a semi-automated brain tumor segmentation method based on regularized non-negative matrix factorization. User-defined seeding points in the pathological regions are combined with a sophisticated seeding method for the normal brain tissues to initialize the NMF algorithm. Piece-wise spatial smoothness as well as sparseness of the NMF tissue abundance maps are encouraged through L1-regularization. Morphological post-processing based on the spatial location of the user-defined seeds is exploited to further remove false positives. The proposed method is applicable to all types of gliomas with any MP-MRI dataset, as the user-defined prior knowledge is patient-specific. We illustrate segmentation performance on an MP-MRI dataset of 21 HGGs combining cMRI, PWI and DWI. To verify robustness against user input variability, each patient’s dataset is repeatedly analyzed with randomly selected seeding points from the pathological subcompartments.


Patient population

Twenty-one patients who were diagnosed with a HGG were enrolled in the study: 1 grade II astrocytoma with focal progression to grade III anaplasia, 2 grade II oligodendrogliomas with focal progression to grade III, 1 oligo-astrocytoma with focal progression to grade III, 6 anaplastic oligo-astrocytomas, 1 anaplastic astrocytoma with focal progression to GBM and 10 GBMs. The Ghent University Hospital local ethics committee allowed a retrospective analysis of the data.

Multi-parametric MRI dataset

The MR examinations were performed on a 3T Siemens Trio Tim scanner (Erlangen, Germany), using a standard 12-channel phased array head coil. All patients underwent an MP-MRI acquisition protocol, consisting of cMRI, PWI and DWI.

cMRI consisted of a 3-dimensional T1-weighted gradient-echo sequence (MPRAGE) before and after contrast administration, with isotropic voxels, and a 3-dimensional T2-weighted inversion recovery sequence (FLAIR) with isotropic voxels. An overview of the image acquisition settings which are defined for all MRI modalities is given in Table 1.

Table 1 Overview of the MR acquisition parameters for cMRI, PWI and DWI

PWI was performed by using a lipid-suppressed, T2*-weighted echo-planar imaging sequence. A series of 90 multi-section acquisitions was acquired at 1 second intervals. The first 10 acquisitions were performed before contrast agent injection to establish a pre-contrast baseline. At the tenth acquisition, a 0.1 mmol/kg body weight bolus of gadobutrol (Gadovist, Bayer) was injected with a power injector (Spectris, Medrad Inc., Indianola, PA) at a rate of 4 ml/s through a 18-gauge intravenous catheter, immediately followed by a 20 ml bolus of sodium chloride solution at 4 ml/s. Relative cerebral blood volume (rCBV) maps were derived from the dynamic signal intensity curves in the DSCoMAN software (Dynamic Susceptibility Contrast MR Analysis, Duke University, Durham, NC). DSCoMAN computes rCBV based on the method proposed by Boxerman et al. [36], compensating for contrast agent leakage due to disruption of the blood-brain barrier. The measured relaxivity change in each voxel is approximated as a linear combination of the whole-brain average relaxivity change in non-enhancing voxels, \(\overline {\Delta R_{2}^{*}(t)}\), and its time integral:

$$ \Delta R_{2}^{*}(t)\approx{K_{1}\overline{\Delta R_{2}^{*}(t)}-K_{2}\int_{0}^{t} \overline{\Delta R_{2}^{*}(\tau)}d\tau} $$

where \(\Delta R_{2}^{*}(t)\) represents the measured relaxivity change in a voxel. The first term reflects the uncontaminated relaxivity change and the second term reflects the effects of leakage. Linear least squares fitting is applied to compute the weighting factors K 1 and K 2 over all the voxels. Corrected relaxivity curves are obtained by only withholding the K 1-term.

Axial diffusion-weighted images were acquired using a fast single-shot gradient-echo echo-planar imaging sequence with diffusion gradient b-values of 0, 500 and 1000 s/mm 2. The b500 and b1000 images were acquired in 3 orthogonal directions. An affine coregistration was applied to account for eddy currents. Apparent Diffusion Coefficient (ADC) maps were derived from the 3 b-values using weighted linear least squares fitting [37]:

$$ \beta = (X^{T}WX)^{-1}X^{T}Wy $$

where β is the diffusion model’s parameter vector, containing the ADC value, X is the design matrix of all diffusion gradients, and y is a vector containing the logarithm of the signal intensities. W is a weighting matrix to take into account the heteroscedasticity, i.e. the fact that the lower signals in y have a higher variance, as a result of the logarithmic transform:

$$ W = diag(exp(2X\beta_{LLS})) $$

where β LLS represents the initial estimate of β obtained using standard linear least squares fitting. The raw b0 images were also added to the input dataset, serving as a T2-weighted reference.

Six MRI features were obtained from the raw acquired data after pre-processing: T1, T1C, FLAIR, rCBV, ADC and b0. All MP-MRI features were coregistered and resampled to the same spatial resolution of 1×1×3 mm 3. cMRI data were skull-stripped and T1C served as a reference for rigid coregistration in SPM8 (Wellcome Trust Centre for Neuroimaging, University College London), using the normalized mutual information criterion [38] and cubic B-spline interpolation for reslicing. Analyses were performed on 10 axial slices located around the tumor centre. Additional intensity-based features were added to the feature set to include localized spatial information. An in-plane local neighborhood of 3×3 and 5×5 voxels was used to calculate average intensity values that were assigned to the central voxel. These spatially averaged intensity values were added for all MP-MRI features. The averaged features were added in accordance with a previous study, in which it was shown that their inclusion significantly improved MRI-based brain tumor segmentation [28]. The same finding was also confirmed for the current study. Each feature’s full range was rescaled linearly to [0–1]. A total of 18 MRI features was finally obtained, making up the rows of the input matrix X.

Segmentation framework

Non-negative matrix factorization

NMF provides a rank-r approximation of a non-negative input matrix X by the product of 2 non-negative factor matrices, W and H:

$$ X\approx{WH} \;\;\text{with}\quad X\in{R^{m\times{n}}_{+}}, W\in{R^{m\times{r}}_{+}} \quad\text{and}\quad H\in{R^{r\times{n}}_{+}} $$

with m being the number of input features and n the number of data points in X. NMF reveals an additive parts-based structure of the input data. It represents each column of X by a weighted sum of the r columns of W. As we are dealing with image intensities, the non-negativity constraint applies naturally. Each column of X corresponds to one voxel’s MP-MRI feature set. Each column of W represents a tissue-specific signature, i.e. an MP-MRI feature vector corresponding to one pure tissue type. As such, each voxel’s MP-MRI feature vector is approximated as a weighted sum of (tissue-specific) source vectors. Each column of H represents the weights of the tissue types for one voxel. One row of H contains the abundances of one particular tissue type over all the voxels, which can be transformed back into the image space to obtain a tissue abundance map. The following objective function is considered for solving the NMF problem:

$$ \min_{W,H}f(W,H)=\min_{W,H}\quad\frac{1}{2}\left(\Vert{X-WH}\Vert^{2}_{F}+\lambda\Vert{(L+I)H}\Vert_{1}\right) $$

The objective function consists of 2 terms. The first term minimizes the difference between the input matrix X and its factorization, WH, based on the Frobenius norm. The second term is an L1-regularization term, which consists of 2 components. The first component, LH, promotes piece-wise smoothness on the tissue abundance maps. L is a sparse n ×n matrix, with each row containing a vectorized Laplacian kernel. As such, each row of L applies a two-dimensional second order spatial derivative to the corresponding voxel. An in-plane neighborhood of 4 voxels is considered for the Laplacian kernel. Due to the relatively low out-of-plane resolution of the MRI data, neighboring voxels in adjacent slices are not considered. The second component applies L1-regularization to H directly, imposing sparseness to the abundance maps. The above NMF formulation [39] is solved using the structured data fusion framework [40] as implemented in Tensorlab [41]. A transformation of variables is used to convert the constrained optimization problem in Eq. 5 into an unconstrained problem, by squaring the entries of the factor matrices. The Gauss-Newton algorithm with dogleg trust region (GNDL) is used to solve the resulting non-linear least-squares problem [42]. At each iteration of the GNDL algorithm, a step is calculated by iteratively solving a linearized version of Eq. 5 using conjugate gradients. A maximum number of 500 iterations and a convergence tolerance of 10−6 for the relative difference of two subsequent values of the objective function were used as stopping criteria for the NMF analyses. The regularization coefficient λ was empirically set to 0.1, after testing a range of λ values and refining the search near λ=0.1.

NMF initialization

As NMF poses a non-convex optimization problem, initialization of the factor matrices (W 0 and H 0) is required. Careful initialization is important, as it might speed up the convergence process and influence the final segmentation result. The pathological tissue sources are initialized based on voxel selection. The user is required to select one or more voxels in each pathological region, i.e. active tumor and, if present, necrosis and edema. For each selected voxel, a candidate source vector is added to W 0 as the averaged feature vector of the selected voxel and its 4 in-plane neighboring voxels. For each pathological tissue class, we then calculate the correlation coefficient between the candidate source vectors as their inner vector product after normalization. The initial pathological sources are obtained after merging candidate source vectors with a correlation coefficient higher than 0.95 in W 0 and replacing them by their average, thereby reducing complexity of the NMF model.

As the normal brain tissue types (i.e. white matter, grey matter, cerebro-spinal fluid and blood vessels) are still highly abundant in the affected brain, an automated seeding procedure is used to obtain their initial source vectors. To cope with the variance within the tissue classes, 2 sources are assigned to each normal tissue type, leading to a total of 8 normal tissue sources. We use the Successive projection algorithm (SPA) [43] to obtain an initial estimate of the normal sources. SPA returns a subset of voxels from the input matrix X with minimal collinearity in feature space. All columns of X are first projected onto the orthogonal subspace of the already initialized pathological sources. A first normal seeding point is selected as the voxel with the highest L2-norm in the orthogonal subspace. All columns are then further projected onto the orthogonal subspace of the added voxel, and the next seeding point is selected as the point with the highest L2-norm among these projected columns. The above procedure is repeated until 8 voxels (for initializing the 8 normal tissue sources) have been selected. Similarly to the pathological sources, the normal sources are obtained as the averaged feature vector of each selected voxel and its 4 in-plane neighboring voxels, assuming these to belong to the same tissue type as the selected voxel. As SPA is known to be sensitive to outliers [44], an additional FCM procedure is applied, using the initial pathological sources and the normal sources obtained from SPA to initialize the cluster centroids. FCM alternatingly updates the cluster centroids and the cluster membership values. As we already obtained proper initialization of the pathological tissue types through user input, the pathological cluster centroids are forced to remain the same throughout the FCM updating procedure. The final FCM centroids make up the columns of W 0. The columns of H 0 are found by applying non-negative least squares fitting to the corresponding column of X and W 0. To reduce computation time, a downsampling factor of 2 was applied to the x- and y-direction of each image slice. After NMF analysis, the obtained tissue abundance maps were upsampled again through linear interpolation. Performance was not affected by this downsampling procedure.

Morphological post-processing

After NMF analysis, each voxel is assigned to the source for which it has the highest abundance value. Voxels that were assigned to different sources belonging to the same pathological tissue class are merged to obtain a preliminary tissue segmentation mask. A morphological post-processing procedure, consisting of 2 steps is then applied to the pathological segmentation masks to further remove false positives, by exploiting the known location of the user-selected voxels. Spatial consistency of the tissue regions is assumed, therefore the segmentation masks are analyzed in terms of (3D-)connected components. For each selected pathological tissue type, only the connected component closest to each user-defined voxel is withheld in the corresponding tissue segmentation mask (Step 1). As the user-defined voxels are not necessarily selected close to the centre of the tissue region, the distance between a selected voxel and a connected component is defined as the minimal distance between the selected voxel and any voxel from the connected component. It is possible that the same connected component is found for several or all user-selected voxels of the same tissue class.

To avoid missing disjoint regions of a tumor component for which no voxels were selected by the user, an additional step is performed which assumes spatial connectivity of the various pathological components (Step 2). For instance, necrotic areas are always adjacent to an active tumor region. Therefore, any connected component of the preliminary necrosis mask which has a common edge with any withheld active tumor component is also included. Similarly, components of the preliminary active tumor mask are also included if they have a common edge with any withheld necrotic component. The same procedure is also applied to edema, by verifying spatial adjacency to the active tumor region(s). Figure 1 illustrates the morphological post-processing procedure for the necrotic tissue mask. One small necrotic component is not withheld after Step 1, as the user didn’t click in this region. However, it was still recovered in Step 2 by verifying spatial connectivity with the active tumor mask.

Fig. 1

Illustration of morphological post-processing after initial semi-automated NMF based segmentation of necrosis a and active tumor c. Step 1: false positives are removed by withholding only the connected components closest to the user-defined seeding points (marked by cursor arrows) for necrosis (b) and for active tumor (d). Step 2: spatial adjacency of the connected components in the preliminary necrosis mask (green) to the withheld active tumor mask (red) is verified in (e). The final necrosis mask is shown in yellow in (f)


Segmentation results of the pathological tissue regions were compared against manual segmentation by an experienced radiologist. The Dice-score was used to quantify the spatial alignment between semi-automated and manual segmentation:

$$ Dice_{\text{{tissue}}}=2\times{\frac{A_{\text{{tissue,NMF}}}\cap{A_{\text{{tissue,man}}}}}{A_{\text{{tissue,NMF}}}+{A_{\text{{tissue,man}}}}}} $$

where A tissue,NMF is the area segmented by NMF and A tissue,man the area manually segmented by the radiologist for the same tissue type. Additionally, the Hausdorff distance was calculated for evaluating the distance between segmentation boundaries. The Hausdorff distance is the maximum distance of all points from one segmentation mask to the corresponding nearest point of the other segmentation mask:

$$ \begin{aligned} Haus_{\text{{tissue}}}&= max\left(\sup_{p\in A_{\text{{tissue,NMF}}}} \inf_{t\in A_{\text{{tissue,man}}}}d(p,t), \sup_{t\in A_{\text{{tissue,man}}}}\right.\\ &\qquad\qquad\left.\inf_{p\in A_{\text{{tissue,NMF}}}}d(t,p)\right) \end{aligned} $$

where sup and inf represent the supremum and infinum, respectively. d(p,t) is a distance metric, for which the Euclidean distance is commonly used. The Hausdorff distance is however susceptible to small outlying subregions in either segmentation masks, as it considers the maximum surface distance. To overcome this limitation, we considered a more robust version of the Hausdorff measure, reporting the 95-percentile instead of the maximum surface distance. In analogy to previous work [14], Dice-scores and Hausdorff distances are reported for active tumor, the tumor core (active tumor + necrosis) and the whole tumor (tumor core + edema). To verify robustness of the semi-automated method to user input variability, the NMF analysis was repeated 20 times per patient with different selection of the seeding points. In each run, an automatic random selection of 3 points was performed in each pathological region, based on the manual segmentation masks.

Besides NMF with spatial regularization and sparseness (\(\mathbf {NMF_{\text {spatial\_sparse}}}\)), additional NMF analyses were performed to assess the added value of the regularization term and the advanced MRI modalities: NMF without regularization (\(\mathbf {NMF_{\text {no\_reg}}}\)), NMF with only spatial regularization (N M F spatial), NMF with spatial regularization and sparseness but without the morphological post-processing step (\(\mathbf {NMF_{\text {no\_postproc}}}\)), and NMF with spatial regularization and sparseness when only considering cMRI data (N M F cMRI). For the N M F cMRI analyses, rCBV and ADC were omitted from the MP-MRI feature set but b0 was withheld, acting as a surrogate for T 2-weighted MRI.


Figure 2 gives an example of the regularized NMF results and tissue segmentation for a GBM patient with a typical ring-enhancing lesion. Some of the MP-MRI input maps are shown on the first row. The NMF abundance maps as well as the final segmentation masks are shown for the pathological tissue types on the second row.

Fig. 2

First row: coregistered MP-MRI maps of a GBM patient, left to right: T1C, FLAIR, rCBV, ADC. Second row, left to right: NMF abundance maps for active tumor, necrosis and edema. The final segmentation masks are shown on the right for active tumor (red), necrosis (yellow) and edema (blue)

Figure 3 gives a comparison of the segmented pathological regions for a GBM patient obtained from the different NMF analyses (\(\mathbf {NMF_{\text {no\_reg}}}\), N M F spatial, N M F cMRI and \(\mathbf {NMF_{\text {spatial\_sparse}}}\)) and manual segmentation. For the active tumor region, spatial overlap with manual segmentation is slightly lower for N M F spatial compared to the other methods. On the other hand, the spatial regularization term did allow N M F spatial to segment the entire necrotic region, whereas the other methods missed a central portion of the necrotic area. Segmentation of edema is inferior for N M F cMRI, where the differentiation from surrounding healthy brain structures was found difficult. \(\mathbf {NMF_{\text {spatial\_sparse}}}\) shows the least false positive regions for necrosis and edema compared to the other methods.

Fig. 3

Comparison of the segmentation results of the pathological tissue regions obtained for a GBM patient using the different NMF methodologies. The left figure of each row shows the obtained segmentation for active tumor (purple), necrosis (red) and edema (yellow). The second to fourth figure of each row show the individual segmentation for active tumor, necrosis and edema, respectively. The top row shows manual segmentation, whereas the other rows show the overlap between the NMF (blue) and manual (green) segmentation result. Segmentation overlap is marked in cyan

Figure 4 shows the dispersion of the Dice-scores for active tumor per patient over 20 runs. Sixteen out of 21 patients have a median Dice-score of at least 60%. Fifteen patients have a lower quartile Dice-score of at least 50%. Seventeen patients have an interquartile range not higher than 15%. Seven patients had at least one run with a Dice-score lower than 30%. A large spread in the Dice-scores is seen for patient 6, who has a non-enhancing anaplastic astrocytoma with focal areas of enhancing GBM. Variability in the results was caused by the difficult differentiation of non-enhancing tumor from edema. Patient 8 suffers from a non-enhancing bifocal GBM, with only a small enhancing area in the smaller lesion. Due to the random selection of voxels in the active tumor region, in some runs voxels were only selected from one of both lesions, resulting in a failure to detect the other lesion. A large dispersion in the Dice-scores is also found for patient 18, exhibiting a heterogeneous and irregularly shaped GBM with varying degrees of enhancement in the active tumor region. Active tumor and necrosis were not well differentiated in some runs, mainly due to ambiguous voxel selection near the pathological tissue boundaries. The low outlier scores for patient 10 are due to the difficult differentiation of non-enhancing tumor from edema. The low outlier score for patient 13 is explained by a suboptimal and unrepresentative random voxel selection in the active tumor region.

Fig. 4

Boxplots showing the dispersion of the Dice-scores for active tumor. Boxplots show quartile ranges of the Dice-scores, ’+’ indicates outliers

Boxplots for the tumor core are shown in Fig. 5. Eighteen out of 21 patients have a median Dice-score higher than 65%. For 17 patients, the lower quartile Dice-score was at least 60%. The interquartile range was not higher than 15% for 17 patients. Seven patients have at least one run with a Dice-score lower than 40%. For 5 out of the 7 patients, these low Dice-scores were found to be outliers (i.e. at a distance of more than 1.5 times the interquartile range from the lower quartile value).

Fig. 5

Boxplots showing the dispersion of the Dice-scores for the tumor core. Boxplots show quartile ranges of the Dice-scores, ’+’ indicates outliers

Figure 6 shows the boxplots for the whole tumor region. A median Dice-score higher than 70% was found for 19 out of 21 patients. Eighteen patients have lower quartile value higher than 60%. Five patients have at least one run with a Dice-score lower than 45%. For 4 out of these 5 patients, such low values were found to be outliers.

Fig. 6

Boxplots showing the dispersion of the Dice-scores for the whole tumor region. Boxplots show quartile ranges of the Dice-scores, ‘+’ indicates outliers

Table 2 gives a comparison of the mean Dice-scores over all patients for the different NMF analyses. Overall, the best performance is obtained with \(\mathbf {NMF_{\text {spatial\_sparse}}}\). When considering only spatial regularization (N M F spatial), segmentation results are mainly worse for the active tumor region, where the Dice-scores are even lower than for \(\mathbf {NMF_{\text {no\_reg}}}\). The spatial smoothing was found to be too severe for several GBMs with a narrow ring-enhancing active tumor compartment. Compared to \(\mathbf {NMF_{\text {no\_reg}}}\), Dice-scores are higher with \(\mathbf {NMF_{\text {spatial\_sparse}}}\) for active tumor and the tumor core, with an increase of 2% in most cases. The Dice-score for the whole tumor region is 0 to 1% higher for \(\mathbf {NMF_{\text {spatial\_sparse}}}\) compared to \(\mathbf {NMF_{\text {no\_reg}}}\). \(\mathbf {NMF_{\text {spatial\_sparse}}}\) does not show any improvement compared to N M F cMRI for the active tumor region. An increase in Dice-score of 1 to 2% is found for the tumor core and an increase of 3 to 4% for the whole tumor region. The lowest Dice-scores are found for \(\mathbf {NMF_{\text {no\_postproc}}}\), with a decrease of 6 to 8% compared to \(\mathbf {NMF_{\text {spatial\_sparse}}}\).

Table 2 Mean Dice-scores [%] for NMF without regularization, with spatial regularization, with spatial and sparse regularization but without morphological post-processing, with spatial and sparse regularization on the cMRI data only, and with spatial and sparse regularization on the full MP-MRI dataset

Table 3 reports the mean Hausdorff distances over all the patients for the different NMF analyses. The lowest Hausdorff distances are found for \(\mathbf {NMF_{\text {spatial\_sparse}}}\), with mean values of 6.1 mm, 7.4 mm and 8.2 mm for active tumor, the tumor core and the whole tumor region, respectively. Comparable but slightly higher Hausdorff distances are found for \(\mathbf {NMF_{\text {no\_reg}}}\) and N M F spatial. Higher Hausdorff distances are found for N M F cMRI, with mean values of 7.4 mm, 9.1 mm and 14.1 mm. Hausdorff distances increase considerably when morphological post-processing is omitted: mean values of 26.6 mm, 29.3 mm and 28.9 mm are found for \(\mathbf {NMF_{\text {no\_postproc}}}\).

Table 3 Mean Hausdorff distances [mm] for NMF without regularization, with spatial regularization, with spatial and sparse regularization but without morphological post-processing, with spatial and sparse regularization on the cMRI data only, and with spatial and sparse regularization on the full MP-MRI dataset


Comparison to other segmentation studies is sometimes hampered by the fact that different segmentation metrics or a different definition of the pathological subregions are being considered. We have quantified our results for active tumor, the tumor core and the whole tumor region using the Dice-score, in accordance with the Multimodal Brain Tumor Segmentation (BRATS) challenge held at the Medical Image Computing and Computer Assisted Intervention (MICCAI) conference. Looking at the Dice-scores reported for HGGs on the BRATS 2012 and 2013 datasets [14], it can be seen that our results are competitive, even when looking at the mean of 25th percentile values (see Table 2). For the whole tumor region, Dice-scores are close to or even higher than 80%, which is in the range of inter-observer variability [14]. We have also applied semi-automated L1-regularized NMF directly to the BRATS 2013 Leaderboard dataset, allowing for a direct comparison with state-of-the-art [45]. It was found that our segmentation framework outperformed all other methods in segmenting the active tumor region, and was also competitive for the tumor core and the whole tumor region. As the BRATS dataset only contains cMRI data, comparison is most appropriate to the N M F cMRI results. The inclusion of additional MRI modalities for brain tumor segmentation has been commonly suggested [2, 46] and initial studies have found improved segmentation results with extended MP-MRI datasets [28, 47]. When including ADC and rCBV into the MRI feature set, we found no advantage for the active tumor region, but Dice-scores increased by 1 to 2% for the tumor core and by 3 to 4% for the whole tumor region. These findings are in accordance with [28], where inclusion of PWI and DWI were found to be mainly advantageous for the tumor core and the whole tumor region. Segmentation performance was also assessed using the Hausdorff distance (see Table 3). As for the Dice-score, we obtained mean Hausdorff distances which are competitive with the best methods reported on the BRATS 2012 and 2013 datasets [14]. Another important consideration of any segmentation algorithm is its computational cost. To reduce computation time, we have applied downsampling to the in-plane dimensions of the imaging data. A downsampling factor of 2 was found not to affect segmentation performance. In the case of isotropic imaging data, where the out-of-plane resolution is as high as the in-plane resolution, downsampling could be applied to all 3 dimensions, with no expected loss in segmentation performance.

Most of the unsupervised brain tumor segmentation algorithms in literature incorporate prior knowledge in the form of basic assumptions or heuristics, which often limits their general applicability to a subset of gliomas or to a particular set of MP-MRI data. Prior knowledge in the form of user input is flexible, allowing to incorporate patient-specific information regarding appearance, location and/or shape of the tumor. Using our semi-automated framework, we were able to tackle the main limitations of many unsupervised classification methods like NMF. Most unsupervised algorithms require proper initialization to come to a valid locally optimal solution. Several studies cope with this by running the segmentation algorithm numerous times with a randomized initialization, then selecting the final solution based on a predefined objective function [24, 28]. We are combining user-defined voxel selection for initializing the pathological sources with a sophisticated initialization of the normal brain tissue sources based on SPA and FCM. As SPA is sensitive to outliers, its output is fed into the FCM algorithm, thereby providing a deterministic initialization for FCM.

Due to the lack of labelled training data, automatic assignment of a tissue label to each segmented region is non-trivial for unsupervised methods. Several studies do not explicitly propose an automated labelling strategy [25, 27], while others assume specific characteristics of the tumor compartments in specific MRI images (such as e.g. contrast enhancement of active tumor) [22, 23]. In our proposed methodology, tissue labelling automatically results from the voxel selection in the pathological regions.

False positive regions were excluded from the pathological tissue masks by exploiting the spatial location of the selected voxels in a morphological post-processing procedure. As the pathological tissue classes are assumed to form spatially consistent regions, only the connected component closest to each selected voxel is withheld in the segmentation masks. Loss of tumor components not selected by the user is avoided by assuming spatial adjacency of the pathological regions. Mean Dice-scores decreased by 6 to 8% when omitting the morphological post-processing step after NMF with spatial regularization and sparseness (Table 2, \(\mathbf {NMF_{\text {no\_postproc}}}\)). Loss in performance was more pronounced when considering the Hausdorff distance, with mean values being an order of magnitude larger for \(\mathbf {NMF_{\text {no\_postproc}}}\) than for the other NMF analyses (Table 3). This indicates that our semi-supervised approach allows for the removal of false positive regions at considerable distance from the actual tumor volume, by exploiting the spatial location of the user-defined seeding points. Similar types of post-processing have been proposed to enhance segmentation results. Cordier et al. only withheld 1 or 2 connected components for each tissue region based on their size, assuming that the largest component(s) correspond to the true tissue region [32]. Menze et al. exploit knowledge about shape and location of the tumor based on an atlas tissue prior to remove false positives [25]. Havaei et al. initialize a k-nearest neighbors (kNN) classifier based on voxel selection by the user [35]. They do not exploit the spatial location of the selected voxels in a post-processing step, but instead they add voxel coordinates to the feature set of the kNN classifier.

Nowadays, a common approach to further improve segmentation results is to model the spatial dependency between the tissue labels of adjacent voxels, typically using a Markov Random Field (MRF) [25, 35] or a Conditional Random Field (CRF) [15, 17]. With MRF, a fixed penalty term is added to the cost function when adjacent voxels have different labels. This penalty term is also depending on the (dis)similarity of the feature vectors when using CRF. To the authors’ knowledge, this is the first study to consider spatially regularized NMF for tumor segmentation. We have added a spatial regularization term with sparseness to the NMF objective function. L1-regularization was used to promote piece-wise smoothness of the abundance maps, allowing for discontinuities at the tissue boundaries. Whereas MRF and CRF are directly applied to the voxel labels, we applied L1-regularization to the tissue abundance maps which represent continuous variables. In some patients, suboptimal segmentation results were found when using spatial regularization only. This is reflected in the lower Dice-scores for N M F spatial in the active tumor region and in the tumor core region (Table 2). The original abundance maps have a relatively low degree of sparseness, such that piece-wise spatial smoothness was not always restricted to the true tissue boundaries. Results were improved by combining spatial regularization with sparseness. This kind of regularization resembles MRF more closely, since MRF implies absolute sparseness of the tissue labels. Over all patients, we found an increase in mean Dice-score by 2% for active tumor and for the tumor core when using spatial regularization with sparseness (Table 2). For the whole tumor region, an increase of 1% was found for the 25th and 75th percentile Dice-scores and no increase was found for the mean value. Similar improvements in segmentation results have been reported in the literature [17, 35].

As imaging characteristics vary considerably across glioma subtypes and grades, supervised classification methods require extensive training datasets with a uniform acquisition protocol. To increase specificity, separate classifiers are usually built for low- and high-grade gliomas [14]. Using training data from different patients to segment a new test case requires careful calibration of image intensities across different patients, which is non-trivial in the presence of pathology. Our method uses labelled data only from the patient under study. Feature vectors averaged over the local neighborhood of the selected voxels are assumed to be approximative prototypes of the pathological tissue regions to be segmented.

One of the main limitations of any semi-automated segmentation algorithm is that they don’t provide reproducible results. Robustness against user input variability is therefore a key aspect. We have assessed user input variability by repeatedly selecting random voxels in each pathological region. Fairly robust segmentation results were found for most patients (see Figs. 4, 5 and 6), but a large spread in the Dice-scores and low outlier values were found in several patients as well. In some cases suboptimal results could be explained by ambiguous voxel selection near the tissue boundaries. Figure 7 illustrates a case of suboptimal voxel selection for patient 1. One of the necrotic voxels was selected at the edge of the manually segmented necrotic region (red voxel in B, other necrotic voxels were selected on other slices). This voxel is located very near the active tumor region and does not show the hypo-intense T1C signal which is characteristic for necrosis. The resulting NMF segmentation is shown in C for active tumor and in D for necrosis. The necrotic voxel selected at the boundary results in an overestimation of the necrotic region and an underestimation of the tumor region. The low outlier Dice-score for patient 1 for active tumor in Fig. 4 corresponds to this voxel selection. Another patient had a bi-focal tumor, and in some runs only one of both lesions was detected as no voxels were selected in the other lesion. Random voxel selection from the manually segmented tumor subregions might not be entirely representative for the user input variability that is to be be expected from trained radiologists. It might be more effective to consider actual user input from various experts in a future study. Valid segmentation results can only be obtained when seeding points are selected in an intelligible way. Spatially distributing the selected voxels is one way of covering intra-tissue heterogeneity. Seeding points should be selected in each lesion in order to properly detect multi-focal tumors. Careful voxel selection is expected to further improve segmentation results and reduce variability.

Fig. 7

Example of a bad segmentation result due to suboptimal voxel selection. A close-up of a GBM lesion on an axial T1C slice (a). The manually segmented necrotic region is shown in yellow (b), the selected necrotic voxel is marked in red. Segmentation of the active tumor region (c) and necrotic region (d) on several slices, blue indicates NMF segmentation, green indicates manual segmentation and cyan indicates overlap


We have presented a semi-automated brain tumor segmentation method, based on NMF with L1-regularization to promote spatial consistency and sparseness of the tissue abundance maps. The semi-automated framework was applied to an MP-MRI dataset consisting of cMRI, PWI and DWI. User-defined voxel selection is applied to initialize pathological sources of the NMF analysis, and to exploit knowledge about the spatial location of the tumor to remove false positives in a post-processing step. In this way, we aimed at incorporating prior knowledge while maintaining general applicability of our method to any type of glioma and to any MP-MRI dataset. Sensitivity to user input variability was explored through repeated analyses with different voxel selection. Robust results were found for most patients, although careful voxel selection is mandatory to avoid sub-optimal segmentation. Based on the reported mean Dice-scores and Hausdorff distances, segmentation results are competitive with state-of-the-art in literature.



Apparent diffusion coefficient


Multimodal Brain Tumor Segmentation challenge


Conventional MRI


Conditional random field


Diffusion-weighted imaging


Fuzzy C-means clustering


Fluid-attenuated inversion recovery


Glioblastoma multiforme


Gauss-Newton with dogleg trust region


High-grade glioma


k-nearest neighbors


Multi-parametric MRI


Markov random field


Non-negative matrix factorization


Perfusion-weighted imaging


Response assessment in neuro-oncology


Relative cerebral blood volume


Successive projection algorithm


T1-imaging with contrast enhancement


Echo time


Inversion time


Repetition time


  1. 1

    Ostrom QT, Gittleman H, Liao P, Rouse C, Chen Y, Dowling J, Wolinsky Y, Kruchko C, Barnholtz-Sloan J. CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2007–2011. Neuro-oncology. 2014; 16(suppl 4):1–63.

    Article  Google Scholar 

  2. 2

    Cha S. Update on brain tumor imaging: from anatomy to physiology. Am J Neuroradiol. 2006; 27(3):475–87.

    CAS  PubMed  Google Scholar 

  3. 3

    Young GS. Advanced MRI of adult brain tumors. Neurol Clin. 2007; 25(4):947–73.

    Article  PubMed  Google Scholar 

  4. 4

    Van Cauter S, De Keyzer F, Sima DM, Sava AC, D’Arco F, Veraart J, Peeters RR, Leemans A, Van Gool S, Wilms G, et al. Integrating diffusion kurtosis imaging, dynamic susceptibility-weighted contrast-enhanced MRI, and short echo time chemical shift imaging for grading gliomas. Neuro-oncology. 2014; 16(7):1010–21.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5

    Cha S, Knopp EA, Johnson G, Wetzel SG, Litt AW, Zagzag D. Intracranial mass lesions: Dynamic contrast-enhanced susceptibility-weighted echo-planar perfusion MR imaging 1. Radiology. 2002; 223(1):11–29.

    Article  PubMed  Google Scholar 

  6. 6

    Bulakbasi N, Kocaoglu M, Farzaliyev A, Tayfun C, Ucoz T, Somuncu I. Assessment of diagnostic accuracy of perfusion MR imaging in primary and metastatic solitary malignant brain tumors. Am J Neuroradiol. 2005; 26(9):2187–199.

    PubMed  Google Scholar 

  7. 7

    Hu LS, Eschbacher JM, Heiserman JE, Dueck AC, Shapiro WR, Liu S, Karis JP, Smith KA, Coons SW, Nakaji P, et al. Reevaluating the imaging definition of tumor progression: perfusion MRI quantifies recurrent glioblastoma tumor fraction, pseudoprogression, and radiation necrosis to predict survival. Neuro-Oncology. 2012; 14(7):919–30.

    Article  PubMed  PubMed Central  Google Scholar 

  8. 8

    Price S, Jena R, Burnet N, Hutchinson P, Dean A, Pena A, Pickard J, Carpenter T, Gillard J. Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: an image-guided biopsy study. Am J Neuroradiol. 2006; 27(9):1969–74.

    CAS  PubMed  Google Scholar 

  9. 9

    Wen PY, Macdonald DR, Reardon DA, Cloughesy TF, Sorensen AG, Galanis E, DeGroot J, Wick W, Gilbert MR, Lassman AB, et al. Updated response assessment criteria for high-grade gliomas: response assessment in neuro-oncology working group. J Clin Oncol. 2010; 28(11):1963–72.

    Article  PubMed  Google Scholar 

  10. 10

    Dhermain F. Radiotherapy of high-grade gliomas: current standards and new concepts, innovations in imaging and radiotherapy, and new therapeutic approaches. Chin J Cancer. 2014; 33(1):16–24.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. 11

    Lacroix M, Abi-Said D, Fourney DR, Gokaslan ZL, Shi W, DeMonte F, Lang FF, McCutcheon IE, Hassenbusch SJ, Holland E, Hess K, Michael C, Miller D, Sawaya R. A multivariate analysis of 416 patients with glioblastoma multiforme: prognosis, extent of resection, and survival. J Neurosurg. 2001; 95(2):190–8.

    CAS  Article  PubMed  Google Scholar 

  12. 12

    Ellingson BM, Cloughesy TF, Lai A, Nghiemphu PL, Mischel PS, Pope WB. Quantitative volumetric analysis of conventional MRI response in recurrent glioblastoma treated with bevacizumab. Neuro-oncol. 2011; 13(4):206.

    Article  Google Scholar 

  13. 13

    Gordillo N, Montseny E, Sobrevilla P. State of the art survey on MRI brain tumor segmentation. Magn Reson Imaging. 2013; 31(8):1426–38.

    Article  PubMed  Google Scholar 

  14. 14

    Menze B, Reyes M, Van Leemput K, et al. The multimodal brain tumor image segmentation benchmark (BRATS). IEEE Trans Med Imaging. 2015; 34(10):1993–2024.

    Article  PubMed  Google Scholar 

  15. 15

    Bauer S, Fejes T, Slotboom J, Wiest R, Nolte LP, Reyes M. Segmentation of brain tumor images based on integrated hierarchical classification and regularization. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer: 2012. p. 10–13.

    Google Scholar 

  16. 16

    Zikic D, Glocker B, Konukoglu E, Criminisi A, Demiralp C, Shotton J, Thomas O, Das T, Jena R, Price S. Decision forests for tissue-specific segmentation of high-grade gliomas in multi-channel MR. In: International Conference on Medical Image Computing and Computer-Assisted Intervention, vol. LCNS 7512. Berlin: Springer: 2012. p. 369–76.

    Google Scholar 

  17. 17

    Lee CH, Wang S, Murtha A, Brown MR, Greiner R. Segmenting brain tumors using pseudo–conditional random fields. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer: 2008. p. 359–66.

    Google Scholar 

  18. 18

    Verma R, Zacharaki EI, Ou Y, Cai H, Chawla S, Lee SK, Melhem ER, Wolf R, Davatzikos C. Multiparametric tissue characterization of brain neoplasms and their recurrence using pattern classification of MR images. Acad Radiol. 2008; 15(8):966–77.

    Article  PubMed  PubMed Central  Google Scholar 

  19. 19

    Jensen TR, Schmainda KM. Computer-aided detection of brain tumor invasion using multiparametric MRI. J Magn Reson Imaging. 2009; 30(3):481–9.

    Article  PubMed  PubMed Central  Google Scholar 

  20. 20

    Li SZ. Markov random field modeling in image analysis: Springer; 2009. ISBN: 1848002793.

  21. 21

    Lafferty J, McCallum A, Pereira F. Conditional random fields: Probabilistic models for segmenting and labeling sequence data. In: Proceedings of the Eighteenth International Conference on Machine Learning, ICML, vol. 1. Williamstown: 2001. p. 282–9.

  22. 22

    Fletcher-Heath LM, Hall LO, Goldgof DB, Murtagh FR. Automatic segmentation of non-enhancing brain tumors in magnetic resonance images. Artif Intell Med. 2001; 21(1):43–63.

    CAS  Article  PubMed  Google Scholar 

  23. 23

    Kazerooni AF, Mohseni M, Rezaei S, Bakhshandehpour G, Rad HS. Multi-parametric (ADC/PWI/T2-w) image fusion approach for accurate semi-automatic segmentation of tumorous regions in glioblastoma multiforme. Magn Reson Mater Phys Biol Med. 2015; 28(1):13–22.

    Article  Google Scholar 

  24. 24

    Juan-Albarracín J, Fuster-Garcia E, Manjón JV, Robles M, Aparici F, Martí-Bonmatí L, García-Gómez JM. Automated glioblastoma segmentation based on a multiparametric structured unsupervised classification. PloS One. 2015; 10(5):0125143.

    Article  Google Scholar 

  25. 25

    Menze BH, Van Leemput K, Lashkari D, Weber MA, Ayache N, Golland P. A generative model for brain tumor segmentation in multi-modal images. In: International Conference on Medical Image Computing and Computer-Assisted Intervention, vol. 13. Berlin: Springer: 2010. p. 151–9.

    Google Scholar 

  26. 26

    Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a hidden markov random field model and the expectation-maximization algorithm. IEEE Trans Med Imaging. 2001; 20(1):45–57.

    CAS  Article  PubMed  Google Scholar 

  27. 27

    Nie J, Xue Z, Liu T, Young GS, Setayesh K, Guo L, Wong ST. Automated brain tumor segmentation using spatial accuracy-weighted hidden markov random field. Comput Med Imaging Graph. 2009; 33(6):431–41.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 28

    Sauwen N, Sima DM, Van Cauter S, Veraart J, Leemans A, Maes F, Himmelreich U, Van Huffel S. Hierarchical non-negative matrix factorization to characterize brain tumor heterogeneity using multi-parametric MRI. NMR Biomed. 2015; 28(12):1599–624.

    CAS  Article  PubMed  Google Scholar 

  29. 29

    Ortega-Martorell S, Lisboa PJ, Vellido A, Simões RV, Pumarola M, Julià-Sapé M, Arús C. Convex non-negative matrix factorization for brain tumor delimitation from MRSI data. PloS One. 2012; 7(10):47824.

    Article  Google Scholar 

  30. 30

    Prastawa M, Bullitt E, Moon N, Van Leemput K, Gerig G. Automatic brain tumor segmentation by subject specific modification of atlas priors 1. Acad Radiol. 2003; 10(12):1341–8.

    Article  PubMed  PubMed Central  Google Scholar 

  31. 31

    Saha BN, Ray N, Greiner R, Murtha A, Zhang H. Quick detection of brain tumors and edemas: A bounding box method using symmetry. Comput Med Imaging Graph. 2012; 36(2):95–107.

    Article  PubMed  Google Scholar 

  32. 32

    Cordier N, Menze B, Delingette H, Ayache N. Patch-based segmentation of brain tissues. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer: 2013. p. 6–17.

    Google Scholar 

  33. 33

    Kwon D, Shinohara RT, Akbari H, Davatzikos C. Combining generative models for multifocal glioma segmentation and registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin: Springer: 2014. p. 763–70.

    Google Scholar 

  34. 34

    Hamamci A, Kucuk N, Karaman K, Engin K, Unal G. Tumor-cut: segmentation of brain tumors on contrast enhanced MR images for radiosurgery applications. IEEE Trans Med Imaging. 2012; 31(3):790–804.

    Article  PubMed  Google Scholar 

  35. 35

    Havaei M, Jodoin PM, Larochelle H. Efficient interactive brain tumor segmentation as within-brain kNN classification. In: ICPR. Stockholm: 2014. p. 556–61.

  36. 36

    Boxerman J, Schmainda K, Weisskoff R. Relative cerebral blood volume maps corrected for contrast agent extravasation significantly correlate with glioma tumor grade, whereas uncorrected maps do not. Am J Neuroradiol. 2006; 27(4):859–67.

    CAS  PubMed  Google Scholar 

  37. 37

    Veraart J, Sijbers J, Sunaert S, Leemans A, Jeurissen B. Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls. NeuroImage. 2013; 81:335–46.

    Article  PubMed  Google Scholar 

  38. 38

    Maes F, Collignon A, Vandermeulen D, Marchal G, Suetens P. Multimodality image registration by maximization of mutual information. IEEE Trans Med Imaging. 1997; 16(2):187–98.

    CAS  Article  PubMed  Google Scholar 

  39. 39

    Wang YX, Zhang YJ. Nonnegative matrix factorization: A comprehensive review. IEEE Trans Knowl Data Eng. 2013; 25(6):1336–53.

    Article  Google Scholar 

  40. 40

    Sorber L, Van Barel M, De Lathauwer L. Structured data fusion. IEEE J Selected Topics Signal Process. 2015; 9(4):586–600.

    Article  Google Scholar 

  41. 41

    Vervliet N, Debals O, Sorber L, Van Barel M, De Lathauwer L. Tensorlab v3.0. 2016. Accessed 24 June 2016.

  42. 42

    Nocedal J, Wright SJ. Numerical Optimization: Springer; 2006. ISBN: 0387400656.

  43. 43

    Araújo MCU, Saldanha TCB, Galvão RKH, Yoneyama T, Chame HC, Visani V. The successive projections algorithm for variable selection in spectroscopic multicomponent analysis. Chemometrics Intell Lab Syst. 2001; 57(2):65–73.

    Article  Google Scholar 

  44. 44

    Gillis N. Successive nonnegative projection algorithm for robust nonnegative blind source separation. SIAM J Imaging Sci. 2014; 7(2):1420–50.

    Article  Google Scholar 

  45. 45

    Sauwen N, Acou M, Sima D, Maes F, Himmelreich U, Achten E, Van Huffel S. A semi-automated segmentation framework for MRI based brain tumor segmentation using regularized nonnegative matrix factorization. In: Proceedings of the 12th International Conference on Signal Image Technology & Internet Based Systems. Naples: IEEE: 2016. p. 88–95.

    Google Scholar 

  46. 46

    Bauer S, Wiest R, Nolte LP, Reyes M. A survey of MRI-based medical image analysis for brain tumor studies. Phys Med Biol. 2013; 58(13):97.

    Article  Google Scholar 

  47. 47

    Di Costanzo A, Scarabino T, Trojsi F, Giannatempo GM, Popolizio T, Catapano D, Bonavita S, Maggialetti N, Tosetti M, Salvolini U, d’Angelo V, Tedeschi G. Multiparametric 3T MR approach to the assessment of cerebral gliomas: tumor extent and malignancy. Neuroradiology. 2006; 48(9):622–31.

    Article  PubMed  Google Scholar 

Download references


Not applicable.


NS, DMS and MA received funding from the Research Foundation Flanders (FWO), grant number G.0869.12N. JV received funding from the Henri Benedictus Fellowship of the Belgian American Educational Foundation. SVH received funding from the government agency for Innovation by Science and Technology (IWT), grant number IM 135005. SVH received funding from the European Research Council under the European Union’s Seventh Framework Programme, grant number 339804. SVH and DMS received funding from the European Research Council under the European Union’s Seventh Framework Programme, grant number 316679. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The datasets analysed during the current study are not publicly available due to them containing information that could compromise research participant privacy/consent. The data are available from the corresponding author on reasonable request.

Authors’ contributions

EA and MA acquired and anonymized the patient imaging data. JV provided software code and advice to pre-process the DWI data. NS and DMS designed the semi-automated segmentation framework. NMF analyses were performed by NS. The study was supervised by SVH and EA. NS, DMS, SVH, FM and UH were major contributors in writing the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Written informed consent was not obtained from all patients for their data to be used in the current study. The need for consent was waived by the local ethics committee of the University Hospital Ghent due to the retrospective nature of the analyses. All patient data were fully anonymized at their native institution prior to access.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information



Corresponding author

Correspondence to Nicolas Sauwen.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sauwen, N., Acou, M., Sima, D. et al. Semi-automated brain tumor segmentation on multi-parametric MRI using regularized non-negative matrix factorization. BMC Med Imaging 17, 29 (2017).

Download citation


  • MRI
  • Segmentation
  • Brain tumors
  • Non-negative matrix factorization
  • Unsupervised classification