The present invention relates to methods and associated devices for analysing magnetic resonance imaging (MRI) data of a tumour to determine the tumour type. The invention is particularly applicable to distinguishing between Glioblastoma multiformes (GBMs) and brain metastases (METs).

GBMs are one of the most common and lethal intracranial tumours. Even with advances in surgical and clinical neuro-oncolog their prognosis remains poor.

METs are another common brain neoplasm in adults. For metastases to occur, cancer cells are released or break off from their primary site, migrate to the central nervous system and develop their own blood supply. In the context of brain metastases, they can lay dormant for various lengths of time before undergoing further growth. If a non-brain primary cancer can be identified early enough, while still localised, a good prognosis may be expected; however, once the tumour has metastasised to the brain, death is inevitable with rare exceptions.

The management of GBM and MET is different. As a result, the accurate differentiation between tumour types using a non-invasive technique may allow early effective treatment whilst minimising the requirement for surgical biopsy, thus affecting prognosis and outcome. In some instances, it is straightforward and uncomplicated to judge whether a tumour is GBM or MET based on the clinical history of systemic cancer or the multiplicity of lesions. Differentiation, however, is always problematic when the lesion is solitary and clinical findings are non-contributory. The heterogeneous structure and morphological appearance of a tumour is related to its cell-type of origin and grade of malignancy; thus, higher classification accuracies have been obtained in differentiating between glioma and meningioma, but much lower classification accuracies (60-80%) have been reported in distinguishing GBM and MET.

GBMs display rapid irregular cellular proliferation and their invasiveness will influence their morphology. It is evident that the shape of a GBM is highly variable. In contrast, METs typically expand more homogeneously. Thus, the shape of METs is expected to be more spherical. This insight has been exploited to distinguish between GBMs and METs. In L. Blanchet, P. W. T. Krooshof, G. J. Postma, A. J. Idema, B. Goraj, A. Heerschap, and L. M. C. Buydens, "Discrimination between Metastasis and Glioblastoma Multiforme Based on Morphometric Analysis of MR Images," American Journal of Neuroradiology, vol. 32, no. 1, pp. 67-73, 2011, MRI data of a tumour was analysed in a slice by slice manner to determine how circular the tumour section in each slice was by determining the proportion of the smallest bounding square around the tumour that is filled by the tumour in each slice. They claimed an accuracy of 93.9% in distinguishing between GBMs and METs for their test data. However, it would be desirable to provide a method capable of improved and/or more reliable accuracy relative to their approach.

It is an object of the invention to provide an improved analysis of MRI data that allows more reliable differentiation between different types of tumour. According to an aspect of the invention, there is provided a method of analysing MRI data of a tumour in order to determine the tumour type, comprising: segmenting the MRI data to identify tumour voxels and non-tumour voxels; processing a plurality of slices of the MRI data in order to obtain, for each slice, values of a plurality of different shape parameters, each shape parameter describing an aspect of the shape of a region containing the tumour voxels within the slice; for each shape parameter, analysing a distribution of values of the shape parameter over the plurality of slices to obtain a value of a shape feature that is characteristic of the tumour shape; and determining the tumour type using the obtained shape feature values.

Thus, a method is provided that uses a plurality of different shape features that are each derived from two-dimensional properties of the shape of a tumour in an MRI slice. The inventors have demonstrated that the use of a plurality of such shape features, particularly when statistically selected according to their effectiveness in distinguishing between the tumour types of interest, allows the determining of tumour type to be carried out more accurately than is possible with prior art techniques that use a single shape feature.

In an embodiment, a combination of shape and spectral features are used to determine tumour type. This approach is particularly effective for distinguishing between GBM and MET tumours. The global morphological appearance of a GBM is expected to have many variances. In contrast, the global shape of METs is expected to be more spheroidal. In addition, both GBM and MET are surrounded by extensive oedema regions on T2 -weighted MR images. The peritumoural oedema around MET is assumed to be increased water content in normal brain tissue, whereas for GBM the oedema region is expected to contain tumour cells that have infiltrated into the tissue. The inventors have found that the combination of 2D shape features, which estimate efficiently the overall exterior shape of the tumour core, and spectral features which estimate efficiently the subtle interior variances of the tumour and oedema regions, is particularly effective in differentiating GBM from solitary MET.

In an embodiment, the determination of the tumour type is based on a plurality of shape features that includes one or more measures dependent on the Solidity of the tumour region, defined as a scalar specifying the proportion of the pixels in the convex hull that are also in the tumour region, computed as -— — , where Area is number of voxels in the tumour region and Convex Area is the number of voxels in the Convex Image. In comparison to a method that determines tumour type using only the Extent of the tumour region, defined as the ratio of the number of voxels in the tumour region to the number of voxels in the rectangular bounding box (which is similar to the approach of Blanchet et al discussed above, differing only in that a rectangular rather than square bounding box is used), better performance is achieved, particularly in distinguishing between GBM and MET. This may occur because Solidity is more sensitive to variations in the heterogeneity of the tumour shapes, which is an important factor in distinguishing between different types of tumour, particularly between GBM (very heterogeneous in shape) and MET (much less heterogeneous in shape). As mentioned above, embodiments can be applied to the problem of distinguishing between METs and GBMs, but they may also be applied effectively to distinguishing between other types of tumour.

The method may also be used to monitor the evolution of a tumour in order to detect a change in the characteristics of the tumour, which may signal the need for a change in treatment or surgery.

According to an alternative aspect of the invention, there is provided a data analysis unit for analysing MRI data of a tumour in order to determine the tumour type, comprising: an input unit configured to receive the MRI data; a data processing unit configured to perform the following steps: segmenting the MRI data to identify tumour voxels and non-tumour voxels; processing a plurality of slices of the MRI data in order to obtain, for each slice, values of a plurality of different shape parameters, each shape parameter describing an aspect of the shape of a region containing the tumour voxels within the slice; for each shape parameter, analysing a distribution of values of the shape parameter over the plurality of slices to obtain a value of a shape feature that is characteristic of the tumour shape; and determining the tumour type using the obtained shape feature values.

According to an alternative aspect of the invention, there is provided an MRI machine, comprising: a data acquisition system for acquiring MRI data of a tumour in a subject; and a data analysis unit for analysing the MRI data of a tumour to determine the tumour type, the data analysis unit comprising: a data processing unit configured to perform the following steps: segmenting the MRI data to identify tumour voxels and non-tumour voxels; processing a plurality of slices of the MRI data in order to obtain, for each slice, values of a plurality of different shape parameters, each shape parameter describing an aspect of the shape of a region containing the tumour voxels within the slice; for each shape parameter, analysing a distribution of values of the shape parameter over the plurality of slices to obtain a value of a shape feature that is characteristic of the tumour shape; and determining the tumour type using the obtained shape feature values.

According to an alternative aspect of the invention, there is provided a computer program for analysing MRI data of a tumour in order to determine the tumour type, the computer program being such that when run on a computer it causes the computer to perform the following steps: segmenting the MRI data to identify tumour voxels and non-tumour voxels; processing a plurality of slices of the MRI data in order to obtain, for each slice, values of a plurality of different shape parameters, each shape parameter describing an aspect of the shape of a region containing the tumour voxels within the slice; for each shape parameter, analysing a distribution of values of the shape parameter over the plurality of slices to obtain a value of a shape feature that is characteristic of the tumour shape; and determining the tumour type using the obtained shape feature values.

Embodiments of the invention will now be described, by way of example only, with reference to the accompanying drawings in which corresponding reference symbols indicate corresponding parts, and in which:

Figure 1 is a flow chart depicting the framework of a method of analysing MRI data of a tumour in order to determine the tumour type;

Figure 2 is a schematic illustration of a tumour region and its rectangular bounding box;

Figure 3 is a schematic illustration of a tumour region comprising a concave region, together with a corresponding convex image and rectangular bounding box;

Figure 4 shows an example tumour region for the purposes of illustrating tumour slice orientation;

Figure 5 shows an ellipse having the same second-moments as the tumour region of Figure 4 superimposed on the tumour region;

Figure 6 depicts the ellipse of Figure 5 in isolation to illustrate the minor and major axes and the slice orientation;

Figure 7 is an artificially generated spectral plot for explanatory purposes;

Figure 8 shows a plot of tumour region spectral features derived from a plurality of GBM cases;

Figure 9 shows a plot of tumour region spectral features derived from a plurality of MET cases;

Figure 10 shows a plot of oedema region spectral features derived from a plurality of GBM cases;

Figure 11 shows a plot of oedema region spectral features derived from a plurality of MET cases;

Figure 12 is a boxplot of 10 exemplary selected shape and spectral features;

Figure 13 depicts ROC curves for a Detailed Example 2;

Figure 14 depicts an MRI machine comprising an MRI data acquisition system comprising a data analysis unit.

According to an embodiment, there is provided a method of analysing magnetic resonance imaging (MRI) data of a tumour in order to determine the tumour type. A framework of the method is illustrated in Figure 1.

In a first step (SI) MRI data is segmented to identify which of the voxels in the MRI data are associated with tumour and which voxels are not. Various approaches for performing this segmentation are applicable to the invention, including approaches using a variety of different MRI modalities or combinations of modalities. The voxels deemed to be tumour voxels may for example be those that are estimated to contain more than a certain proportion of cancerous cells. Similarly, the voxels deemed to be non-tumour voxels may be those that are estimated to contain less than a certain proportion of cancerous cells. Other metrics for distinguishing between "tumour voxels" and "non-tumour voxels" may be used. What is important is that the segmentation between the tumour voxels and non-tumour voxels enables a volume of tissue to be identified that has a surface that broadly represents a boundary between tissue that is (mostly) cancerous and tissue that is (mostly) normal (or at least mostly non-cancerous). In an embodiment, the segmentation is performed manually, for example by inspecting a visual representation of the MRI data and indicating manually, based on clinical experience, where a boundary between tumour and non-tumour voxels is thought to lie.

In an alternative embodiment, segmentation is performed using an automated process (i.e. carried out by a computer), which is time-efficient for tumour delineation whilst minimising inter-observer error. The segmentation may be applied to anatomic MRI alone. Alternatively, the segmentation may operate on quantified diffusion data (DTI) directly.

In an embodiment, the so-called "D-SEG" method is used, which involves generation of a visual display of tumour isotropic and anisotropic diffusion characteristics using p:q maps with minimal observer inputs. The D-SEG method can find the 3D contours, which delimit tissue boundaries between tumour and normal brain by segmentation in the p:q space. In an embodiment, D-SEG applies a -means clustering algorithm, which iteratively segments p:q space into k non-overlapping clusters. The -means clustering defines a prototype in terms of centroids and is applied to objects in a continuous n-dimensional space, i.e., p:q space. The first step of -means clustering is to define k initial centroids, the number of which is specified a priori according to the number of clusters desired. In MRI image segmentation, this decision is necessarily made based on functional and anatomical considerations. In an embodiment, k is set to be 16, which has been determined to identify particularly effectively the range of potential tissue compartments present within a brain affected by a tumour, e.g., normal appearing white matter, normal appearing grey matter, cerebro-spinal fluid spaces, solid tumour, regional tumour necrosis, tumour-associated cystic regions, peri-lesional oedema, peri-lesional tumour infiltration and distant oedema. Secondly, space is separated into k subsets based on percentiles of the p and q data present within the scans. Initial cluster centroids are determined as the median coordinate in p:q space for each cluster. Next, the distance is calculated from each voxel to each cluster centroid in p:q space. Each voxel is then assigned to its nearest cluster in p:q space and cluster centroids, i.e., medians, are then recalculated based on the new data within these clusters. This procedure is iteratively repeated until no point changes clusters, the centroids remain the same or a defined iterative limit is

where ( ;, q ) represents the voxel in the p:q space, and (m _p , m _q ) denotes its nearest cluster median. The degree of scatter of each cluster or measure of distance of data points to its corresponding centroid is measured as a sum of the squared error. Eventually, the flood-fill method is employed by putting seed voxels inside the tumour region to semi-automatically generate regions of interest (ROIs) slice by slice.

Binary masks of the tumour core region, i.e. the tumour shape region, can then be obtained.

Additionally, binary masks of the oedema regions can also be obtained. In a further step (S2) the segmented MRI data is processed in a slice by slice manner in order to calculate values of a plurality of shape features. The shape features are derived from shape parameters, as follows.

For each slice, values of a plurality of different shape parameters are calculated. Each of the shape parameters describes an aspect of the shape of a region containing the tumour voxels within the slice. For example, one or more of the following shape parameters may be used: Area, Rectangular Bounding Box Area, Convex Area, Eccentricity, Equivalent Diameter, Perimeter, Extent, Solidity, Major Axis Length, Minor Axis Length, and Slice Orientation. These shape parameters are defined as follows:

a. Area'— Scalar; the actual number of voxels in the tumour region or "region of interest" (ROI). The ROI may be defined by D-SEG segmentation, for example, and it may also be referred to as a binary mask.

b. 'Rectangular Bounding Box Area'— The smallest rectangle containing the ROI. This is illustrated in Figure 2. Here, an ROI 2 is defined as the region within the perimeter 6 defined by the path ABCDA. The rectangular bounding box 4 is shown and intersects each of the points A, B, C and D, in this example.

c. 'Convex Area'— Scalar that specifies the number of voxels in the 'Convex Image'. The 'Convex Image' is a binary image that specifies the convex hull, with all voxels within the hull filled in (i.e., set to on). The image is the size of the bounding box 4 of the region. This is illustrated in Figure 3. Here, a different ROI 2 is defined as the region within the perimeter 6 defined by the path ABCDEA. Thus, a region 8, which was part of the ROI 2 in the example of Figure 2 is no longer part of the ROI 2 in the example of Figure 3. Nevertheless, it can be seen that the Rectangular Bounding Box Area is the same as in the example of Figure 2 because the rectangular bounding box 4 is identical. ROI 2 is an example of what would be referred to as a concave region. The convex image of the ROI 2 is the region within the perimeter defined by the path ABCDA. Thus, the Convex Area of the ROI 2 in Figure 3 is the same as the Convex Area of the ROI 2 in Figure 2 despite the ROIs being different. d. 'Eccentricity'— Scalar that specifies the eccentricity of the ellipse that has the same second- moments as the ROI. This is illustrated in Figures 4-6. Figure 4 shows an ROI 2. Figure 5 shows an ellipse 10 having the same second-moments as the ROI 2 overlaid on the ROI 2. The ellipse has a major axis 12 and a minor axis 14. Figure 6 shows the ellipse 10 in isolation (for clarity), and shows the ellipse's foci 16. The eccentricity is the ratio of the distance between the foci 16 of the ellipse 10 and its major axis length. The value is between 0 and 1. (0 and 1 are degenerate cases; an ellipse whose eccentricity is 0 is actually a circle, while an ellipse whose eccentricity is 1 is a line segment.) e. 'Equivalent Diameter'— Scalar that specifies the diameter of a circle with the same area as the ROI.

Here: Area'— Scalar; the actual number of voxels in the region.

f. 'Perimeter'— Scalar; the distance around the boundary of the ROI. g. 'Extent'— Scalar that specifies the ratio of voxels in the region to voxels in the total bounding box. Computed as the Area divided by the area of the rectangular bounding box (refer to Fig.2 for the definition of the rectangular bounding box).

h. 'Solidity'— Scalar specifying the proportion of the voxels in the convex hull that are also in the ROI

(refer to Figure 3 for the definition of the convex hull). Computed as— — .

° ConvexArea

i. 'Major Axis Length'— Scalar specifying the length (in voxels) of the major axis of the ellipse that has the same normalized second central moments as the region (see Figures 4-6).

j. 'Minor Axis Length'— Scalar; the length (in voxels) of the minor axis of the ellipse that has the same normalized second central moments as the region (see Figures 4-6).

k. 'Slice Orientations'— Scalar; the angle (in degrees ranging from -90 to 90 degrees) between the x- axis and the major axis of the ellipse that has the same second-moments as the region. This is illustrated in Figures 4-6. In Figure 6, the foci 16 of the ellipse 10 are shown, together with the angle

18 between the x-axis 20 and the major axis 12.

For each shape parameter, a distribution of values of the shape parameter over the plurality of slices is analysed in order to obtain a value of a shape feature that is characteristic of the tumour shape. The shape feature may comprise the maximum, the mean, the standard deviation, the median, or one or more intermediate quantiles such as the 0.25 or 0.75 quantiles, of one of the shape parameters described above, for example. For the Slice Orientations, the orientation in degrees is arbitrary without a fixed reference system, so only the variance (measured by standard deviation) is used in this case.

Therefore, in this example, a total of 61 shape features are available.

In an embodiment, a feature selection process is performed to determine a subset of the 61 shape features that would be most effective for distinguishing between the different types of tumour of interest. Step S2 is then performed using only the shape features selected in the feature selection process.

In an embodiment, the feature selection process comprises analysing a database of MRI data consisting of data obtained from a plurality of different patients having known types of tumour. Values of a plurality of different shape features are obtained from this data for each of a plurality of different tumours of known types. For example, values of a plurality of different shape features are obtained from each of a plurality of different tumours that are all of a known first type and values of the same plurality of different shape features are obtained from each of a plurality of different tumours that are all of a known second type (different from the first type).

For each shape feature, a test is then applied to determine the statistical significance of differences between a set of values of the shape feature obtained from the plurality of tumours of the first type and a set of values of the shape feature obtained from the plurality of tumours of the second type. The test may comprise a Student's T-Test for example. A subset of shape features may then be selected based on the determined statistical significance. For example, this selection may be carried out on the basis of selecting only those shape features where the determined significance is greater than a predetermined threshold value.

The feature selection process may further take into account correlation between different shape features. In an embodiment, for example, pairs of shape features are identified in the selected subset of shape features that are highly correlated with each other. Using both members of any such pair would be inefficient as the information content of each member is very similar. Efficiency can therefore be improved by removing one of the shape features from one or more of the identified pairs in order to form a smaller subset of features for use in step S2.

In an embodiment, at least one shape feature is selected that is obtained by analysing a distribution of values of the Solidity, the Solidity being defined as a scalar that specifies the proportion of the voxels in the convex hull of the tumour region that are also in the tumour region. The Solidity is particularly sensitive to the degree of outline variability of the tumour region and is particularly effective for distinguishing between GBM and MET tumours for example.

In an optional further step (S3), spectral information of either or both of the tumour core region and the tumour with oedema region can be derived. The spectral information comprises one or more spectral features. Each spectral feature may comprise information about the frequency of occurrence of one of the different types of segment derived using the segmentation of step S 1. For example, in the case where the segmentation is performed using the D-SEG method described above, 16 different types of segment may be used, so there could be 16 different types of spectral feature in this example. More generally, such spectral information can be useful whenever the MRI data is segmented into three or more different types of segment.

In an embodiment, each spectral feature is defined as the ratio between the number voxels that are of a particular segment type (e.g. which fall within the segment in p.q space where the D-SEG method is used) and the total number of all voxels in the ROI.

If there are 16 different segments and both tumour and oedema spectra are obtained, this provides a further 32 features that can be used to distinguish between the different tumour types.

A feature selection process that is analogous to the feature selection process for selecting shape features only, discussed above, may be used to select the most effective spectral features to use. In an embodiment, the feature selection process selects the most effective combination of shape and spectral features. The feature selection process may comprise obtaining values of a plurality of spectral features from MRI data obtained from each of a plurality of tumours of a first type and a plurality of tumours of a second type. For each spectral feature, a test may be applied to determine the statistical significance of differences between a set of values of the spectral feature obtained from the plurality of tumours of the first type and a set of values of the spectral feature obtained from the plurality of tumours of the second type. A subset of spectral features may then be selected according to the determined statistical significance. The feature selection process may further comprise identifying pairs of spectral features in the selected subset of spectral features that are highly correlated with each other, and removing one of the spectral features from one or more of the identified pairs in order to form a smaller subset of features.

In step S4, values of the shape features calculated in step S2 are used to determine the tumour type of the MRI data. Detailed Example 1 discussed below illustrates an embodiment of this type. Alternatively, in step S4 a combination of values of the shape features calculated in step S2 and values of the spectral features calculated in step S3 are used to determine the tumour type of the MRI data. Detailed Example 2 discussed below illustrates an embodiment of this type.

The determination of the tumour type in step S4 may be carried out using a supervised learning strategy, for example. A training set for the supervised learning may consist of the MRI data for the tumours of known type used in above-described feature selection steps, for example.

As shown in the Detailed Example 1 below, the inventors have discovered that in the case where shape features only are used for the determination of the tumour type, and in particular for determining whether a tumour is a GBM or MET tumour, the following shape features are particularly effective: a quantile over the slices between 0 and 0.5 of the Extent, a maximum over the slices of the Eccentricity, a maximum over the slices of the Perimeter, a standard deviation over the slices of the Extent, a maximum over the slices of the Solidity. In an embodiment, two or more of these particular shape features are used. In an embodiment, all of these shape features are used. In an embodiment, all of these shape features and no other shape features are used.

As shown in the Detailed Example 2 below, the inventors have discovered that in the case where a combination of shape features and spectral features are used for the determination of the tumour type, and in particular for determining whether a tumour is a GBM or MET tumour, the following shape features are particularly effective: maximum over the slices of the Rectangular Bounding Box Area, a median over the slices of the Extent, a median over the slices of the Solidity, a maximum over the slices of the Solidity. In an embodiment, two or more of these particular shape features are used. In an embodiment, all of these shape features are used. In an embodiment, all of these shape features and no other shape features are used.

Detailed Examples 1 and 2 are now discussed in turn.

Feature Selection, Classification and Performance - Detailed Example 1

In an embodiment, a feature selection process is applied to the shape features only.

A specific example of such an embodiment is described below.

The data acquisition was as follows. MR data were acquired on a GE Signa Horizon 1.5T MR system (GE Healthcare, Milwaukee, WI, USA) equipped with 22mT/m gradients and using a quadrature head coil. DTI data sets were obtained with a diffusion weighted spin echo echo-planar-imaging sequence, that is, one acquisition without diffusion sensitization, i.e., b = 0 sec/mm ² , and one acquisition with diffusion weighting, i.e., b = 1000 sec/mm ² in 12 gradient directions. Coverage of whole brain was with 50 contiguous slices with two interleaved series of four repeats (2.8mm thick slices with 2.8mm gaps, TR/TE = 7000/80ms, acquisition matrix = 96x96, and FOV = 24cm), as described in T. R. Barrick and C. A. Clark, "Singularities in diffusion tensor fields and their relevance in white matter fiber tractography," Neurolmage, vol. 22, no. 2, pp. 481-491, 2004.

The patient subjects were as follows. In agreement with the local regional ethics committee, 48 patients, who were recruited in two blocks between 2005-2006 and 2008-2010, were retrospectively entered into this study. All the patients were not subject to any treatment before the radiologic examinations, and they had no prior history of surgery, chemotherapy, or radiation therapy applied to the intracranial structures. In addition, histopathologic analysis of the resected tissue confirmed the diagnosis of GBM in 30 patients. For the 18 patients with MET, histopathologic analysis shows that 6 cases originated from the lung, 2 cases from breast cancer, 2 from squamous cell, 2 from bowel adenocarcinoma, 1 from melanoma, 1 from renal cancer, 1 from neuroendocrine, 1 from prostate cancer and 2 from unknown source.

In this embodiment, a Student's T-Test and a correlation test is performed in order to determine which combination of the shape features would be most effective for distinguishing between different types of tumour. Features can be ranked according to the /^-values of the statistic hypothesis test, with smaller p- values representing features which more effectively distinguish between different tumour types.

In an example, the inventors applied this process to the case of distinguishing between GBM and MET. By setting the significant level to 0.01, 33 significant features were selected. In a subsequent step, highly correlated features, characterised by absolute values of the Pearson's linear correlation coefficient larger than 0.7, were removed. After these two steps of feature selection procedure, eventually there are five features remaining including the 0.25 quantile of the Extent (qEx ), the maximum Eccentricity (maxEcc), the maximum Perimeter (maxPer.), the standard deviation of the Extent (stdExt), and the maximum Solidity (maxSol.).

In an embodiment, a supervised learning strategy may be used to distinguish between different tumour types (e.g. between GBM and MET) using the selected features. In the present example, the inventors used quadratic linear discriminant analysis (QDA) (P. A. Lachenbruch and M. Goldstein, "Discriminant analysis," Biometrics, vol. 35, no. 1, pp. 69-85, 1979), Naive Bayes (NB) (G. John and P. Langley, "Estimating continuous distributions in Bayesian classifiers," Proceedings of the Eleventh conference on Uncertainty in artificial intelligence, pp. 338-345, 1995), k-Nearest Neighbour algorithm (kNN) (J. Friedman, J. Bentley, and R. Finkel, "An algorithm for finding best matches in logarithmic expected time," ACM Transactions on Mathematical Software (TOMS), vol. 1549, no. July, pp. 209-226, 1977), nonlinear Support Vector Machine (SVM) (N. Cristianini and J. Shawe-Taylor, An Introduction to Support Vector Machines and Other Kernel-based Learning Methods, First. Cambridge University Press, 2000) and Neural Networks (NNWs) (M. T. Hagan, H. B. Demuth, and M. H. Beale., Neural network design. PWS Pub. Co., 1996) to distinguish between GBM and MET.

QDA is a generative classification model based on the multivariate Gaussian. It constructs a statistical framework to calculate the likelihood of target occurrence in future trials given the occurrences in prior trials. In addition, QDA incorporates different covariance matrices for each class that can result in a superior description of individual populations by adding more parameters.

Another popular probabilistic classifier is NB, which estimates the probability an object belongs to a specific class if the observed feature values of an object and the prior probabilities of classes are available. The predicted class is that with the highest estimated probability. The model is "naive" as the features are not assumed to be independent, even conditional on the class label. However, even if the assumption of the NB is not true, it often results in accurate classification because the model is simple and hence it is relatively immune to overfitting.

kNN is a non-parametric classifier separating objects based on nearest training examples in a feature space. It is a memory- or instance-based learning method using a distance or similarity function, e.g., Euclidean distance was applied. kNN classifies an object by estimating an empirical fraction, i.e., a majority vote, of its neighbours, with the object being assigned to the most frequent class among its ^-nearest neighbours. Furthermore, we find the optimal number of k equal to 10 by using cross validation.

SVM is a powerful technique for binary classification. The predictions of SVM only depend on a subset of training data, i.e., support vectors, and it combines with the kernel trick to find the best hyperplane with the largest margin between two classes. In the present example, a Gaussian Radial Basis Function kernel was used with a default scaling-factor, sigma, of 1 to map the feature vectors into a nonlinear feature space. Then an optimal hyperplane is constructed to separate all the data points to either GBM or MET class.

NNWs consist of simple elements operating in parallel and these elements are inspired by biological nervous systems. A series of logistic regression models (LRMs) stacked on top of each other with the final layer being another LRM can form a feed-forward neural network, a.k.a., multi-layer perceptron (MLP) classifier. Thus, the logistic function is used as the transfer function. The architecture of the MLP network is defined using the cross-validation method that the number of hidden layers is 9 according to the best cross- validation accuracy.

In the present example, the classification performance (i.e. the ability to distinguish between GBM and MET in this example) was evaluated by (i) Leave-One-Out (LOO) cross validation, which is an unbiased predictor and is capable of creating sufficient training data for studies with small sample size (B. Efron, "Estimating the error rate of a prediction rule: improvement on cross-validation," Journal of the American Statistical Association, vol. 78, no. 382, pp. 316-331, 1983) (ii) the area under the receiver operating characteristic (ROC) curve (see C. Metz, "ROC methodology in radiologic imaging," Investigative radiology, vol. 21, pp. 720-733, 1986 and H. Zweig, "Receiver-Operating Clinical Medicine ( ROC ) Plots: A Fundamental Evaluation Tool in," vol. 39, no. 4, pp. 561-577, 1993.) and (iii) the classification accuracy, sensitivity and specificity at the operating point on the ROC curve using various classifiers mentioned above.

Table 1 tabulates the /^-values of the two-sample Student T-Test and correlation coefficients of the two steps feature selection.

Table 1 : Feature selection by the two-sample Student T-Test and correlation analysis. Here we tabulate the p-values of the finally selected 5 features and correlation coefficients (Corr coef:

Correlation Coefficients).

uF.xt maxF.cc ma Per SJclF.xt maxSol

T-Test p-values 3.87E-11 4.21E-07 3.06E-06 1.38E-05 1.58E-05

(C oi r coef. (GBM) 1 - - - -

-0.43 1 - - -

0.13 -0.13 1 - -

-0.60 0.46 -0.24 1 -

0.08 0.08 -0.07 0.55 1

Corr coef. (MET) 1 - - - -

-0.41 1 - - -

-0.27 0.30 1 - -

-0.06 0.14 -0.17 1 -

0.57 -0.13 -0.36 -0.04 1

As aforementioned, in this example 5 features were chosen with smallest /^-values and lowest correlation in between. The obtained /^-values of these 5 features are much smaller (ranging from 3.87E-1 1 to 1.58E-05) than the significant level, i.e., 1E-04, that was set. After the second step of the feature selection, the highly correlated features were filtered out. The 5 non-redundant features remaining have a maximum correlation coefficient of -0.6, which is between the qExt and stdExt, which were derived from the same quantity named Extent. In addition, the minimum correlation is -0.04 between stdExt and maxSol.

The numerical results of cross-validation (LOO), i.e., accuracy, sensitivity, specificity, and AUC are summarized in Table 2.

Table 2: Quantitative results of various classification methods on 5 selected shape features (CVAccu., CVSens., CVSpec, and CVAUC: Leave-One-Out Cross Validation Accuracy, Sensitivity, Specificity, and Area Under ROC Curve. Results in parentheses produced using the "Extent" feature only).

Classifiers CVAccu. CVSens. CVSpec. CVAUC

QDA 87.5 (87.5) 83.33 (88.89) 90 (86.67) 0.969 (0.961)

NB 93.75 (89.58) 94.44 (94.44) 93.33 (86.67) 0.985 (0.954) kNN 93.75 (91.67) 94.44 (100) 93.33 (86.67) 0.983 (0.967)

SVM 91.67 (89.58) 94.44 (100) 90 (83.33) 0.959 (0.965)

NNW 97.92 (91.67) 100 (94.74) 96.88 (89.66) 0.996 (0.886)

The results of various classifiers using the 5 selected features were compared. In addition, a comparison was also made with the results of using only one feature ("Extent"; values given in parentheses). Firstly, using 5 selected features, QDA obtained 88% accuracy with a promising AUC (0.969), but the sensitivity is much lower than the specificity. NB and kNN have got the same accuracy (94%) and the same sensitivity and specificity. However, NB had slightly higher AUC compared to kNN. SVM achieved comparable results as NB and kNN with lower specificity and accuracy. NNW outperformed other classifier and offered 98% accuracy with 100% sensitivity and 97% specificity, and also resulted in an AUC of 0.996. Secondly, the result of using only the Extent feature obtained an accuracy up to 92% (here qExt was used because it is highly correlated with both mean and median Extent). Compared to the method of using 5 selected shape features, the classification using Extent only could achieve higher sensitivity except using NNW, but it always resulted in inferior specificity. The AUC of this example method is superior to the method of using Extent only with an exception of using SVM (Table 2).

Feature Selection, Classification and Performance - Detailed Example 2

In an alternative embodiment, a feature selection process is applied to both the shape features and the spectral features.

A specific example of such an embodiment is described below.

The data acquisition was as follows. MR data were acquired on a GE Signa Horizon 1.5T MR system (GE Healthcare, Milwaukee, WI, USA) equipped with 22 mT/m gradients and using a quadrature head coil. DTI data sets were obtained with a diffusion weighted spin echo echo-planar-imaging sequence, that is, one acquisition without diffusion sensitization, i.e., b = 0 sec/mm ² , and one acquisition with diffusion weighted, i.e., b = 1000 sec/mm ² in 12 gradient directions. In addition, the coverage of whole brain contains 50 contiguous slices with two interleaved series of four repeats (2.8mm thick slices with 2.8mm gaps, TR/TE = 7000/80ms, acquisition matrix = 96 X 96, and FOV = 24 cm).

The patient subjects were as follows. In agreement with the local regional ethics committee, 47 patients, who were recruited in two blocks between 2005-2006 and 2008-2010, were retrospectively entered into this study. All the patients were not subjected to any treatments before the radiologic examinations, and they had no prior history of surgery, chemotherapy, or radiation therapy applied to the intracranial structures. In addition, histopathologic analysis of the resected tissue confirmed the diagnosis of GBM in 30 patients. For the 17 patients with MET, histopathologic analysis shows that there are 6 cases originated from the lung, 2 cases from breast cancer, 1 from squamous cell, 2 from bowel adenocarcinoma, 1 from melanoma, 1 from renal cancer, 1 from neuroendocrine, 1 from prostate cancer and 2 from unknown source.

In this embodiment, a filter-based feature selection method is performed, which relies on a relevance index criterion. In an embodiment, the relevant index criterion is estimated by mutual information maximisation (MIM) (see Brown G, Pocock A, Zhao M, Lujan M. Conditional likelihood maximisation: a unifying framework for information theoretic feature selection. The Journal of Machine Learning Research. 2012;13:27-66).

In the specific example discussed below, a feature selection process using a conditional MIM criterion proposed by Fleuret (Fleuret F. Fast binary feature selection with conditional mutual information. The Journal of Machine Learning Research. 2004;5: 1531-1555) was used, which resulted in the selection of 10 features (4 shape features and 6 spectral features). These selected features included maximum Bounding Box (f3), median of the Extent (f4), median of the Solidity (f5), maximum Solidity (flO), three spectral features of the tumour core area (segments 10 (f6), 12 (f7) and 14 (f8)) and three spectral features of the tumour with oedema region (segments 5 (fl), 9 (f2) and 14 (f9)). The segment numbers refer to the segments as shown in Figures 8-11, discussed further below.

In the example, 1 -sample Kolmogorov-Smirnow test (KST) was used to assess normality of the selected features. The KST returns a test decision for the null hypothesis that the selected features come from a standard normal distribution, against the alternative that features are not from such a distribution. In addition, Pearson's χ ² goodness-of-fit (GoF) test was applied to verify if these features are random samples from a normal distribution with mean and variance estimated directly from the features. To compare differences of the features in two tumour groups, a two-sample t-test if the feature follows a normal distribution or a nonparametric Mann- Whitney U-test (MWUT) otherwise was used.

Statistical analysis was performed on commercial statistical package (Statistical Toolbox of Matlab 2013a, Mathworks Inc., Natick, MA, USA), and p < 0.05 was considered statistically significant.

In this example, an ensemble learning based adaptive boosting (AdaBoost) method was used for classification (see Freund Y, Schapire R. A decision-theoretic generalization of on-line learning and an application to boosting. Computational learning theory. 1995) (in this example, for distinguishing in particular between GBM and MET). The AdaBoost method yields a powerful ensemble by combining several weak classifiers, e.g., Adaboost with Decision Trees (Dietterich T. An experimental comparison of three methods for constructing ensembles of decision trees: Bagging, boosting, and randomization. Machine learning. 2000: 139-157), Neural Networks (see Dietterich T. Ensemble methods in machine learning. In: Multiple classifier systems.; 2000: 1-15 and Ratsch G, Onoda T, Miiller K. Soft margins for AdaBoost. Machine learning. 2001:287-320) or Support Vector Machine (SVM) (Li X, Wang L, Sung E. AdaBoost with SVM-based component classifiers. Engineering Applications of Artificial Intelligence. 2008;21(5):785-795). In the boosting method, each weak classifier is successively incorporated and trained on a subset of the training data and tries to reduce the bias of the combined model. In particular, according to previous research, AdaBoost with SVM offers superior performance over its counterparts for unbalanced dataset. Since the dataset of the present example, which contains GBM and MET, is unbalanced, for AdaBoost with SVM method was chosen (see Piatt J. Fast Training of Support Vector Machines using Sequential Minimal Optimization. In: Advances in Kernel Methods - Support Vector Learning. MIT Press; 1998:41-64 and Keerthi SS, Shevade SK, Bhattacharyya C, Murthy KRK. Improvements to Piatt's SMO Algorithm for SVM Classifier Design. Neural Computation. 2001; 13(3):637-649).

The classification performance was evaluated by (i) stratified Leave-One-Out (LOO) cross validation, which is an unbiased predictor and is capable of creating sufficient training data for studies with small sample size, (ii) the area under the receiver operating characteristic (ROC) curve (AUC) and (iii) the classification accuracy at the operating point on the ROC curve using AdaBoost with SVM classifier aforementioned. In addition, several other metrics were introduced to compare the performance of our method with those of others. For each tumour type, the following were computed: the true positive (TP) rate, false positive (FP) rate, precision, F-measure, which is the harmonic mean of precision and recall (recall is equivalent to the TP rate), and the Cohen's κ coefficient that measures the agreement of prediction with the true class— 1.0 signifies complete agreement.

Each of the voxels in the whole brain image can be assigned to one of the 16 segments (from p-q segmentation) and spectral features can then be derived from the statistics of the 16 segments, for example to show the extent to which the tumour is build up from voxels of each of the different tissue types corresponding to the different segments. Further details of the segmentation method may be found in Jones T, Bell B, Barrick T. A novel whole-brain DTI segmentation technique for brain tumour delineation and diagnosis. In: Proceedings of the international society for magnetic resonance in medicine (ISMRM).Vol 20.; 2012: 188, for example. A spectral feature may be illustrated by a spectral plot. Figure 7 is an artificially generated spectral plot (not generated from real data), for the purposes of explanation. The plot illustrates an imaginary situation where a tumour area (segmented for example by further semi-automated flood-fill method from the whole brain segmentation) is determined to contain in total 500 voxels, with 100 of these belonging to segment 10, 100 to segment 12, 250 to segment 16 and 50 to other segments. The vertical axis of the plot is a unitless percentage representing the median of the segmented region composition (as a %) and the horizontal axis is the segment number (between 1 and 16). It can be seen that segment 16 is the highest point at 50% and segments 10 and 12 are at 25%. In this example only a single tumour case is considered. However, if we have many cases to consider, it is possible to obtain median, upper and lower quartiles of the composition (as in Figures 8-11 discussed below).

Figures 8-11 show spectral plots that illustrate the relative contribution of each of the 16 spectral features for tumour and oedema regions and for GBM and MET cases (i.e. the proportion of voxels that belong to each of the 16 different segments in each case). Figure 8 shows tumour region spectral features for GBM cases, Figure 9 shows tumour region spectral features for MET cases, Figure 10 shows oedema region spectral features for GBM cases and Figure 11 shows oedema region spectral features for MET cases. In each case, a solid line is provided that connects the median values, vertical bars above and below the median values represent the bounds of the upper and lower quartiles, and the whiskers and dots represent outliers. Figure 12 is a boxplot of the selected 10 shape and spectral features selected using the conditional MIM method. In this example, the 10 features comprise three spectral features of the tumour with oedema region fl (corresponding to segment 5), f2 (corresponding to segment 9) and f9 (corresponding to segment 14), three spectral features of the tumour core area f6 (corresponding to segment 10), £7 (corresponding to segment 12), f8 (corresponding to segment 14), and 4 shape features f3 (maximum Bounding Box), f4 (median of the Extent), f5 (median of the Solidity), and f 10 (maximum of the Solidity).

KST rejected the null hypothesis for all the features that they do not follow standard normal distribution; however, GoF test showed that three features, i.e., fl, f4 and £5, are normally distributed. Furthermore, t-test showed that there are significant differences of f4 and £5 between GBM and MET (p < 1 X 10 ^~6 ), but f 1 (p = 0.0508) is not significant different between two types of tumours. Additionally, MWUT for other non-normality features showed that flO has significant differences (p < 1 X 10 ^~5 ) between two groups, but not for others (Figure 12). ROC curves are shown in Figure 13 and the AUC for both GBM and MET are 0.927 with an overall accuracy of 95.74% using stratified LOO cross validation (Table 3). TP rate for the GBM is 100% and 88.2% for the MET, and FP rate averages to 0.075 for the two groups. Precision and F-measure are greater than 0.95 averaged, and the overall κ statistic is 0.905.

Table 3: Numerical results of the classification between GBM and MET using AdaBoost with SVM.

Class TP Rate FP Rate Precision F-Measure AUC Accuracy K statistic

GBM 1 0.118 0.938 0.968 0.927

MET 0.882 0 1 0.938 0.927 95.74% 0.905

Weighted Avg. 0.957 0.075 0.960 0.957 0.927

The methods described above could in principle be used to characterise the shape of any lesion and/or to monitor changes in the shape with time due for example to disease progression or treatment response.

At least one of the steps of the method can be performed using a suitably programmed computer comprising hardware such as CPU and RAM with which the skilled person would be very familiar.

As shown in Figure 14, in an embodiment, a data analysis unit 60 is provided that may optionally form part of an MRI machine 62 (or be separate from it). The data analysis unit 60 may comprise an input unit 64 for receiving MRI data, directly or indirectly, from an MRI data acquisition system 66. The data analysis unit 60 may further comprise a data processing unit 68 (comprising standard computer hardware for example) that is configured to carry out the steps of the method of analysing MRI data of a tumour in order to determine the tumour type discussed above. A computer program may be provided for example (e.g. via a computer medium storing the computer program or via a network connection) that comprises instructions which, when run on a computer, for example the data processing unit 68, cause the computer to carry out the steps of the method of analysing MRI data of a tumour in order to determine the tumour type discussed above.