WO2015040434A1 - Analysing mri data to determine tumour type - Google Patents

Abstract

There is provided a method of analysing MRI data of a tumour in order to determine the tumour type. In one arrangement, the method comprises segmenting the MRI data to identify tumour voxels and non-tumour voxels. A tumour surface is identified representing a boundary between the tumour voxels and the non- tumour voxels. The identified tumour surface is processed to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape. A tumour type is determined based on the one or more calculated shape parameters.

Description

ANALYSING MRI DATA TO DETERMINE TUMOUR TYPE

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 (Tsuchiya et al. 2005). Differentiation, however, is always problematic when the lesion is solitary and clinical findings are non-contributory (Law et al. 2002). The heterogeneous structure and morphological appearance of a tumour is related to its cell-type of origin and grade of malignancy; thus, higher classification accuracy have been obtained in differentiating between glioma and meningioma, but much lower classification accuracies (60-80%) have been reported in distinguishing GBM and MET (Devos et al. 2004; Opstad et al. 2004; Georgiadis et al. 2008).

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 (Frieboes et al. 2006). Thus, the shape of METs is expected to be more spherical. This insight has been exploited to distinguish between GBMs and METs. In (Blanchet et al. 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%, but their model is oversimplified and may not be optimally accurate nor optimally reliable. Furthermore, the analysis in terms of data slices has shortcomings in terms of computation efficiency. Moreover, their approach discards top and bottom slices, which can contain useful morphological information.

It is an object of the invention to provide an improved analysis of MRI data that achieves greater accuracy, greater reliability and/or which is more computationally efficient. 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; identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels; processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape;

determining the tumour type based on the one or more calculated shape parameters.

Thus, a method is provided that uses one or more shape parameters that each represents a three- dimensional characteristic of the tumour shape. It has been recognized previously that the shape of different types of tumour are systematically different. For example, it is known that a MET will tend to be more spherical than a GBM. However, previous methods for numerically analysing MRI data of a tumour have focussed on deriving shaped parameters that are each derived from a two-dimensional slice of the tumour. Such a shape parameter does not therefore represent a three-dimensional characteristic of the tumour shape. For example, previous analyses have derived shape parameters that each represents the extent to which a two-dimensional slice fills the smallest square bounding box. A greater space filling is seen on average for spherical tumours in comparison to non-spherical tumours, and this can be used determine tumour type. A plurality of slices may be obtained, which provides an improved analysis of the tumour, but this approach is computationally inefficient in comparison to the present invention (where individual metrics are each able to represent a three-dimensional characteristic of the surface). Furthermore, extreme slices are discarded in such methods, so not the entire tumour is sampled. This involves loss of potentially valuable information.

Furthermore, the analysis of slices is dependent on the direction of slicing and can lead to further loss of valuable information about the three-dimensional shape. Consider for example a tumour shaped like a cucumber. If the cucumber is sliced along the axis of the cucumber all of the slices may have a uniform cross-section that is roughly circular. Deriving shape parameters that measure the extent to which each slice fills the smallest square bounding box would yield a result that is the same as if the tumour was spherical. It is clear that such methods are prone to losing valuable information about the tumour shape. The present invention does not use slices and is not prone to the same shortcomings.

According to an embodiment, the method further comprises identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein: one or more of the shape parameters is/are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches.

This embodiment is based on the inventors' realisation that the nature of the local three-dimensional shape of the tumour surface, and/or variations in this local three-dimensional shape over the tumour surface, can also be characteristic of the tumour type. This information is encoded efficiently by deriving shape parameters that are based on distributions of local three-dimensional shape characteristics of surface patches over the tumour surface. As shown below, use of such shape parameters achieves more accurate results than approaches which do not consider local three-dimensional shape information.

Various measures of the three-dimensional shape ("shape parameters") that are based on distributions of local shape characteristics are available. As will be demonstrated below, the inventors have recognised that the following are particularly effective: average Shape Index (avgSI); standard deviation of the Shape Index ( stdSI); mean Shape Index (mSI); standard deviation of the Curvedness (stdCv); mean Curvedness (mCvN); mean curvature L2-norm (mCN). In an embodiment one or more of these are used.

In an embodiment, the one or more shape parameters that are each representative of a three- dimensional characteristic of the surface shape comprise(s) one or more shape parameters that is/are not derived based on a distribution of local shape characteristics over a plurality of surface patches. For example, a global measure such as the volume ratio of the tumour with respect to the minimum cuboid bounding box, vCR, and/or the volume ratio of the tumour with respect to the minimum convex hull, vCxR, may be used. This approach can improve the accuracy of the determination (classification) process, particularly when used in combination with one or more shape parameters that are derived based on a distribution of local shape characteristics over a plurality of surface patches.

In an embodiment, the one or more shape parameters that are derived based on a distribution of local shape characteristics over a plurality of surface patches comprise avgSI or stdSI but not both, mCN or mCvN but not both, and the one or more shape parameters that are not derived based on a distribution of local shape characteristics over a plurality of surface patches comprise vCxR or vCR but not both. Using a limited number of measures reduces the computational load associated with the analysis and these particular combinations of choices yield particularly accurate results.

In an embodiment, the combination of vCxR, mCvN and stdSI is used, which provides particularly accurate results.

The method may be applied to the problem of distinguishing between METs and GBMs and/or 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; identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels; processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape; determining the tumour type based on the one or more calculated shape parameters.

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; identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels; and processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape; determining the tumour type based on the one or more calculated shape parameters.

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; identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels; processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape; determining the tumour type based on the one or more calculated shape parameters.

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 illustrating the framework of a method according to an embodiment of the invention;

Figure 2 depicts a triangular mesh generated from a binary mask volume that has been segmented manually for an MET test case;

Figure 3 depicts a triangular mesh generated from a binary mask volume that has been segmented manually for a GBM test case;

Figure 4 depicts a triangular mesh generated from a binary mask volume that has been segmented using a D-SEG algorithm for an MET test case;

Figure 5 depicts a triangular mesh generated from a binary mask volume that has been segmented using the D-SEG algorithm for a GBM test case;

Figure 6 illustrates the definition of the principle curvatures of Koenderink et al. 1992. The figure illustrates a surface patch expressed by L = f(x, y) = x² - y² , which is a saddle shape. By rotating the normal plane that is perpendicular to the tangent plane and is pivoted with respect to the normal vector, it is possible to obtain curves by cutting the surface patch using this normal plane. In addition, the curvatures of these resulting curves are called normal curvatures at the origin p . In general, the maximum and minimum of these normal curvatures, i.e., κ_γ and κ₂ at a given point, e.g., ?, on a surface are called the principal curvatures;

Figure 7 illustrates a polar coordinate system of the principal curvature (K_U K₂) plane, in which the Shape Index and the Curvedness are defined. The antipodal points A and C denote opposite shapes, and points A and B are related by a reflection that represent features of 90 degrees-rotated versions of each other;

Figure 8 illustrates a map of the Shape Index overlaid on the triangular mesh in the case where segmentation is performed manually for an MET test case;

Figure 9 illustrates a map of the Shape Index overlaid on the triangular mesh in the case where the segmentation is performed manually for a GBM test case;

Figure 10 illustrates a map of the Shape Index overlaid on the triangular mesh in the case where segmentation is performed using the D-SEG algorithm for an MET test case;

Figure 11 illustrates a map of the Shape Index overlaid on the triangular mesh in the case where the segmentation is performed using the D-SEG algorithm for a GBM test case;

Figure 12 illustrates a map of the Curvedness overlaid on the triangular mesh in the case where segmentation is performed manually for an MET test case;

Figure 13 illustrates a map of the Curvedness overlaid on the triangular mesh in the case where the segmentation is performed manually for a GBM test case;

Figure 14 illustrates a map of the Curvedness overlaid on the triangular mesh in the case where segmentation is performed using the D-SEG algorithm for an MET test case;

Figure 15 illustrates a map of the Curvedness overlaid on the triangular mesh in the case where the segmentation is performed using the D-SEG algorithm for a GBM test case;

Figure 16 is a boxplot comparing the ability of different shape parameters to distinguish between MET and GBM;

Figure 17 is a table showing the ranking of different shape parameters according to their ability to distinguish between MET and GBM derived using T-Test, KS-Test and MWW-Test;

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

The data discussed below corresponds to a specific example of the invention in which MRI data were acquired with a GE Signa Horizon 1.5T MRI system (GE Healthcare, Milwaukee, WI, USA). DTI data sets were obtained using a diffusion-weighted spin echo echo-planar-imaging sequence (Barrick and Clark 2004). In agreement with the local regional ethics committee, 37 patients, recruited in two blocks between 2005-2006 and 2008-2010, were retrospectively entered into the study. 15 GBM cases and 8 MET cases were randomly chosen for the purposes of demonstrating the effectiveness of the example.

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 for an example 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.

Figures 2 and 3, for example, show meshes that have been generated using such a manual process.

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 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 to identify 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 infiltratration and distant oedema. Secondly, p:q space is separated into k subsets based on either quantiles or 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 reached that is expressed as

min/e_{i,2 16} (Pi ^{_ n}p;)² + - ^;)²,

where ( j, c j) represents the voxel in the 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.

Figures 4 and 5, for example, show meshes generated based on segmentation using the D-SEG algorithm.

In a further step (S2) the segmented MRI data is used to identify the tumour surface representing the boundary between the tumour voxels and the non-tumour voxels.

The standard Marching Cubes (MC) method may be used to identify the tumour surface. The overall process may be referred to as "surface extraction" and "mesh generation". The MC method is a well-known volume visualisation method originally investigated by Lorensen and Cline (Lorensen and Cline 1987). It takes as input a regular scalar volumetric data set and creates triangle models of isosurfaces of a scalar function given by samples over a cuberille grid (Lewiner and Lopes 2003; Newman and Yi 2006). Essentially, the MC method is considered as a surface rendering approach for medical image data (Lorensen and Cline 1987; Gong and Zhao 2010), e.g., MRI data, which can be described in a discrete 3D regular data field in which the value of each voxel is defined as follows, fijk = f (xi_> y_r Zk) _> S- t. : i = 1,2, ... , N_x;j = 1,2, ... , N_y; k = 1,2, ... , N_Z.

The MC method subdivides the data space into a series of small cubes and constructs a facetised isosurface by processing the data set in a sequential scanline, e.g., cube by cube. The cubes are defined by lattice of the volume and each lattice point is a corner vertex (V ) of a cube. In addition, the isosurface of interest is given by

S_f ≡ {{x, y, z) \f{x, y, z) = a], in which a is a pre-defined isovalue (Nielson 2003; Gong and Zhao 2010). During the scanline processing, each cube vertex Vi that has a value equal to or above the isovalue is marked and all other vertices are left unmarked (Newman and Yi 2006). The isosurface intersects each cube edge Ej terminated by one marked vertex and one unmarked vertex and a cube is defined as active when it contains an intersected edge. Therefore, there are 2 , i.e., 256, possible marking scenarios for a cube, which has eight vertices that can be either marked or unmarked (Chernyaev 1995). Then each cube-marking scenario can encode a cube based isosurface intersection pattern (Newman and Yi 2006). Due to the reflective and rotational symmetry property of the cube, the standard MC method simplifies the encoding scheme, which results in only 15 unique marking scenarios. In addition, the MC method can store the facetisation information, e.g., the vertices of the triangulated shapes and the edges they intersected, with respect to the 15 or a higher number of basic intersection topologies for resolving internal ambiguity in a look-up table, which can be built offline prior to the processing of marching cubes (Chernyaev 1995; Lewiner and Lopes 2003; Nielson 2003; Newman and Yi 2006).

In an embodiment, according to the facetisation pattern in the look-up table of the intersection topologies, the MC method creates triangular facets, which denote the portion of the isosurface intersected by each cube. A linear interpolation method may then be applied to calculate the isosurface intersection with the cube edge to allow the isosurface-edge intersection locations to be estimated with subvertex accuracy (Newman and Yi 2006). Eventually, given a binary volume created from segmentation, a separating surface can be extracted as an isosurface corresponding to the isovalue a, and all the obtained triangular facets can produce a triangular "mesh". Examples of such meshes are shown in Figures 2-5.

In a further step (S3) the identified tumour surface is processed to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape (e.g. the whole surface shape). In an embodiment, one or more of the shape parameters are based on a distribution of local three-dimensional shape characteristics of a plurality of local surface patches. For example, a plurality of surface patches representing different portions of the tumour surface may be identified and one or more characteristics of the shape of each of these surface patches may be derived. The identified surface patches may together cover all of the tumour surface or a portion consisting of less than the entire tumour surface.

In an embodiment, the local three-dimensional shape characteristics are based on the Shape Index (SI) or Curvedness (Cv) of the surface patches. In (Koenderink and Van Doom 1992), Koenderink and van Doom derived the local shape and curvedness according to the established curvature measures in differential geometry theory. For differential geometry of 3D shapes, there is no global coordinate system. Therefore, assuming that a surface S in 9¾³, where 5 is a local surface patch, is given by the graph of a smooth function L = (x, y), the tangent plane of this surface S is defined as the L = 0 plane, and n = (0,0,1) is a unit normal vector to the surface S, with respect to the origin p of the surface (as shown in Figure 6). In general, at every point p on a differentiable surface in 3 -dimensional Euclidean space, we can derive its unit normal vector. In addition, a normal plane at p is the one that contains the normal; therefore, it also contains a unique direction tangent to the surface and cuts the surface in a plane curve. This curve generally has different curvatures for various normal planes at p. The two principal curvatures κ₁ and κ₂, which measure the maximum and minimum bending of the surface, are the smallest and largest values of these curvatures. Mathematically, for second-order geometry, the 2D isophote landscape form the second order structure matrix, or the local Hessian matrix (Ter Haar Romeny 2011), is as follows

dydx dy²

in which L_xy = L_yx and this Hessian matrix is symmetric. In addition, the eigenvalues of this Hessian are found by solv

which are the two principal curvatures (Ter Haar Romeny 2011). Furthermore, two other local quantities, i.e., the Gaussian curvature Q and the mean curvature K " , are widely used in standard differential geometry. The Gaussian curvature and the mean curvature are the geometrical mean and arithmetic mean of the principal curvatures, respectively (Lin 2003). Gaussian curvature is calculated by the product of the principal curvatures that is equal to the determinant of the Hessian matrix that is

Q = K_tK₂ = det(H) = L_XXL yy ^>xy

and the mean curvature, which is equal to the trace of the Hessian matrix, is simplified as

K = = trace (H)

2 ' 2

Due to all local approximations for which the ratio of the principal curvatures is equal are of the same shape (Koenderink and Van Doom 1992), an intuitive notion of shape can be apparently defined in a polar coordinate system (κ_1# κ₂) (Figure 7). The direction in this polar coordinate system encodes the shape, and the distance from the origin encodes the size. Subsequently, the Shape Index can be defined using the two principal curvatures (Ter Haar Romeny 2011) as follows:

in which κ₁≥ κ₂ , and the Curvedness is defined as the length of the shape vector

Furthermore, the Curvedness is inversely proportional to the size of the shape, and the SI defines a continuous distribution of surface types ranging from a spherical cup-like shape to spherical cap-like shape. Thus, it can be shown that while the SI is invariant by homothety, the Curvedness is not, and therefore, shape information and size can be easily decoupled (Koenderink and Van Doom 1992).

Thus, the Shape Index and Curvedness can each be considered as an example of a local three- dimensional shape characteristic of a surface patch. A plurality of different shape parameters can be derived based on different distributions and/or different Shape Indexes and/or Curvednesses.

For example, the average Shape Index (Awate et al. 2008) may be used as a shape parameter, which is referred to herein as "avgSI" and defined as follows,

The avgSI may be calculated with respect to the i-th vertex of the extracted surface S.

As a further example, the standard deviation of the Shape Index may be used as a shape parameter, which is referred to herein as "stdSI" and defined as follows, stdSI - \^{ηζ ~ m}

where m is the number of surface patch vertices being considered (e.g. where the surface patches are triangular, each patch will have 3 vertices and m represents the total number of vertices in the mesh). As a further example, the mean Shape Index may be used as a shape parameter, which is referred to herein as "mSI" and defined as follows,

mSI = ¾ m .

Alternatively or additionally, it is possible to derive features from the Curvedness that can be used as alternative or additional shape parameters. For example, the average Curvedness may be used as a shape parameter, which is referred to herein as "avgCv" and defined as follows,

avgCv = [/._6S ( (i)d5 ]^a5.

As a further example, the standard deviation of the Curvedness may be used as a shape parameter, which is referred to herein "stdCv" and defined as follows, stdCv =

As a further example, the mean Curvedness may be used as a shape parameter, which is referred to herein as "mCvN" and defined as follows,

mCvN = .

As a further example, the intrinsic curvature index (Awate et al. 2008) may be used as a shape parameter, which is referred to herein as "iCv" and defined as follows,

iCv = f_iesg⁺(0ds ,

where Q⁺ = max (Q, 0).

As a further example, the mean curvature L2-norm may be used as a shape parameter, which is referred to herein as "mCN" and defined as follows,

m N = _ieS H² (i ds ⁵ .

Shape parameters that are not based on a distribution of local shape characteristics of patches may be used additionally or alternatively. For example, a standard deviation of a Normalized Radial Length (stdCD), a volume ratio of the tumour relative to a minimal cuboid bounding box (vCR), and/or a volume ratio of the tumour relative to a minimum convex hull (vCxR) may be derived. These are described in further detail below.

Defining B as the volume of interest (VOI) of the tumour identified in the segmentation step S I, and Sf as the tumour surface (e.g. an "isosurface" of B with discrete vertex value ½ following the terminology discussed above in the context of the MC method); the Normalized Radial Length (NRL) is defined by the function d_n(V) as

°-^{n K}-^l) - max[d(i)] '

in which d(i) = VO(0 - ¾)² + C O ^{~ Y}oY + 0(0 - Z₀)² , 1≤i≤N. (X₀, Y₀, Z₀) and (x(i),y(0_< ^z(0) are the coordinates of the mass centre of B and the i-th vertex on the isosurface Sf, respectively. N is the number of vertices ½ . The distance between the centroid of the VOI is normalised by the maximum value of the radial length max[d(i)] .

The standard devia ion of the NRL (stdCD) can be expressed as follows: stdCD = - d )²,

where d_n(i) is the mean value of d_n(i), which can be interpreted as the radius of a perfect sphere volume that has the same volume of the irregular shape of the tumour.

It is also possible to derive the volume ratio of the tumour with respect to the minimum cuboid bounding box B_s (vCR) and/or the minimum convex hull B_c (vCxR). These features can be derived as follows:

r_t Volume(B)

vCR = — - ,

Volume(B_s)

vCxR = ^Ξ ^ _t

Volume(Bc)

wherein the function Volume(-) furnishes the number of voxels of B, B_s and B_c. The bounding box B_s may be defined according to the axis-aligned bounding box model that the extreme corners delimit the body diagonal of the box, resulting in B_s being given by all (Xs_> s_> ^zs) coordinates satisfying x_mi_n < x_s <

Xmax? y-min — — Vinax ^min — ¾ — ^max-

The convex hull B_c of a vertex set ½ is the smallest convex region containing the vertices. Each edge of the boundary of B_c is a line segment that can be expressed as an implicit linear equation, and the half space containing the hull is given by an inequality ax_c + by_c + cz_c + d < 0 in 3D, where a, b, c, d are floating point numbers. The region inside B_c is defined by the collection of all (x_c, y_c, z_c) coordinates satisfying these inequalities.

In an embodiment, a combination of two or more local three-dimensional shape characteristics is determined for each of the plurality of surface patches (e.g. for each patch, a value of the Shape Index, SI, and a value of the Curvedness, Cv, may be obtained). In such an embodiment, at least one of the shape parameters may be derived based on a distribution of the combination of two or more local three- dimensional shape characteristics over the plurality of surface patches (e.g. based on how many of the patches comprise each possible combination of values of SI and Cv). For example, a density distribution or histogram representing the distribution of the combination of two or more three-dimensional shape characteristics over the plurality of surface patches may be obtained. The density distribution or histogram may be considered as representing a "spectrum" of the types of shape found over the tumour surface. As described below, the spectrum may be represented in a simpler form by dividing the space defining the possible combinations of two or more local three-dimensional shape characteristics into a plurality of regions. In an embodiment, a 2D histogram of the distribution of a combination of two shape features over the whole tumour may be determined. In this 2D example, the combination of two or more local three- dimensional shape characteristics comprises a combination of exactly two local three-dimensional shape characteristics. In such an embodiment, any pair of local three-dimensional shape characteristics can be used. It has been found that the Shape Index (SI) and Curvedness (Cv) are particularly effective. Each point in the 2D histogram will indicate the number of points (e.g. surface patches) on the surface of each tumour with a particular shape (i.e. a particular combination of the two shape characteristics being considered, for example a particular value of SI and a particular value of Cv for the surface patch in question). A method such as k-means clustering can then be used to divide the 2D space into separate regions. For each tumour the normalised number of points within each region then provides a shape parameter that describes the tumour. In this case the shape parameter may take the form of a vector containing a scalar value for each of the regions. For example, if there are six regions, the vector representing the shape parameter may comprise six values. In a simple illustrative example, if 10% of the surface patches fall into the first region, 20% into the second region, none into the third region, 40% into the fourth region, 5% into the fifth region, and 25% into the sixth region, the vector will be as follows: (10%, 20%, 0, 40%, 5%, 25%). This shape parameter (expressed as a vector) may be referred to as a spectrum for the patient. By comparing the shape parameter for the patient with distributions of the spectra obtained previously for different types of tumour, it is possible to determine (with relatively high accuracy) which type of tumour the patient has. The inventors have found for example that the distribution of spectra for GBM cases is significantly different to the distribution of spectra for MET cases (the spectra for different GBM cases tend on average to be more similar to each other than to spectra for MET cases). Applying this approach to a specific example, the inventors found that using two clusters from a k-means analysis with N=6 produced a high cross-validation accuracy of 96%, an AUC of 0.98 and balanced error rate of 4%.

In step (S4), the tumour is classified (i.e. the tumour type is determined/estimated) using at least one of the calculated shape parameters. The potential effectiveness of the method can be appreciated by inspecting Figures 8-11 and Figures

12-15.

Figures 8 and 10 show the Shape Index mapped onto the tumour surface (respectively for manual and D-SEG segmentation) for an MET test case. The relatively uniform shading shows that there is relatively little variation in the Shape Index over the surface. Figures 9 and 11, on the other hand, show the Shape Index mapped onto the tumour surface (respectively for manual and D-SEG segmentation) for a GBM test case. Here, a much greater degree of variation in the Shape Index over the surface of the tumour can be seen. Thus, a measure of the three-dimensional shape of the tumour derived from a distribution of the Shape Index values over the tumour surface can be effective to distinguish between the two types of tumour (MET and GBM).

Figures 12 and 14 show the Curvedness mapped onto the tumour surface (respectively for manual and D-SEG segmentation) for an MET test case. Here, relatively slowly changing shading indicates a correspondingly smooth and gradual variation in Curvedness over the tumour surface. Figures 13 and 15, on the other hand, show the Curvedness mapped onto the tumour surface (respectively for manual and D-SEG segmentation) for a GBM test case. Here, it can be seen qualitatively that the shading varies in a different manner to the shading in the MET case, indicating that the distribution of the Curvedness over the surface is fundamentally different between the two types of tumour and that a shape parameter derived from the distribution of the Curvedness over the tumour surface can be effective to distinguish between the two types of tumour.

A feature selection (FS) process can be performed to determine which shape parameters (which may also be referred to as "shape features") of the three-dimensional shape will be most effective for distinguishing between different types of tumour. The merits of FS are manifold: first, selecting features can avoid overfitting and therefore improve the performance of the models; second, reducing the number of features used can provide more cost-effective models; third, knowledge of which features it is useful to retain can reveal a deeper insight into the underlying processes (Saeys, Inza, and Larranaga 2007).

FS methods designed with various evaluation criteria generally fall into three categories: the filter model, the wrapper model and the hybrid model (Liu and Yu 2005).

Firstly, the filter model attempts to select features based on some auxiliary criteria, e.g., feature correlations, to select a subset of features that avoid redundancy (Levner 2005). The filter methods decouple the FS process, independently of the performance evaluation component of the classification, by selecting subsets of variables as a pre-processing step that is separated from the choice of the predictor; therefore, these methods are tractable and cost-efficient. However, these methods may ultimately select irrelevant features as a result (Levner 2005). Moreover, when there are some features that are not useful by themselves but they can provide valuable information by combining with others, the filter methods show a lack of such information. Essentially, filter methods are designed for a specific type of feature. Instead of selecting features by invoking the classifier using the wrapper methods, filter methods rank individual features in advance. In the present example the inventors used statistical paired difference tests to evaluate if each feature (shape parameter) is differentially expressed between the two classes (i.e. GBM and MET). Paired Student's T-Test, the paired Kolmogorov-Smirnov Test (KS-Test), and the paired Mann-Whitney-Wilcoxon Test (MWW-Test) algorithms, which are commonly used statistics, were used.

Secondly, wrapper methods try to utilise the learning machine of interest as a black box to score subsets of variables based on the predictive power of a predictor (Haury, Gestraud, and Vert 2011). As testing all combinations of variables is computationally impossible, these methods usually perform a greedy search in the subspace of the features. In this work the inventors employed heuristic sequential forward and backward FS approaches, which are known wrapper methods.

In the case of the paired Student's T-test, the null hypothesis is that shape parameter values (the values of avgSI, stdSI, etc.) from both classes have an identical mean. In the case of paired KS-Test, the null hypothesis is that the cumulative distribution of shape parameter values for class GBM is the same as the cumulative distribution of the shape parameter values for class MET. In addition, using the MWW-Test, it is possible to decide whether the population distributions are identical, which is the null hypothesis, without the normality assumptions. All the three tests determine if the observed differences are statistically significant and return a score representing the probability that the null hypothesis is true (Levner 2005). Therefore, shape parameters can be ranked using one of these statistics according to the significance score of each shape parameter. However, instead of extracting the information from only one test with a cut-off threshold, in the present work the inventors ordered the shape parameters according to the /^-values given by the three tests (see table in Figure 17). Figure 16 shows a boxplot of selected shape parameters versus normalised values of the shape parameters. The shaded boxes 50 show the interquartile range for GBM and the shaded boxes 52 show the interquartile range for MET. Broadly speaking, where these boxes 50 and 52 are well separated from each other, e.g. for vCxR or mCvN, the associated shape parameters are likely to be effective for distinguishing GBM from MET. In contrast, where the boxes 50 and 52 are closer together or overlapping, e.g. avgCv or iCv, the associated shape parameters may be less useful for distinguishing GBM from MET. Circled points are outliers (examples of such points are labelled 54).

From both the boxplot (Figure 16) and the /^-values in Figure 17, it can be seen that the best shape parameters, which are the most separable between two classes, are vCR, vCxR, mCN, mCvN, stdSI and avgSI, and are highlighted by the shaded boxes 56 in Figures 17. In addition, features such as iCv and avgCv, which have large /^-values derived from all the three tests are less effective in discriminating between GBM and MET.

Furthermore, Pearson product-moment correlation coefficient (r) is a well-known and effective measure of similarity between two random variables. Using this measure the inventors have observed that vCR and vCxR are highly correlated (r = 0.86, p < 0.001 for all the GBM cases and r = 0.86, p < 0.001 for all the MET cases), mCN and mCvN are also highly correlated (r = l, p < 0.001 for all the GBM cases and the same for all the MET cases), and large correlation coefficients are obtained for stdSI and avgSI (r = 0.95, p < 0.001 for all the GBM cases and r = 0.97, p < 0.001 for all the MET cases). Therefore, the method can be implemented effectively if only one shape parameter is chosen from each of the three groups {vCr, vCxR}, {mCN, mCvN} and {stdSI, avgSI}, e.g. vCxR, mCvN and stdSI.

The inventors applied three pattern classification methods for comparison: LDA with Fisher's linear discriminant analysis (FLDA) (Fisher 1936), k-nearest neighbour (k-NN) (Friedman, Bentley, and Finkel 1977), and nonlinear support vector machines (SVMs) (Suykens and Vandewalle 1999; Cristianini and Shawe-Taylor 2000). First, the FLDA method is searching for a linear combination of the variables that has a maximum between-class difference relative to its within-class standard deviation. Essentially, the objective of the FLDA method is finding a transformation that minimises the overlapping of the transformed distributions because usually there is no complete separation provided by the transformation (Lachenbruch and Goldstein 1979). Second, the k-NN is a method for classification based on closest training examples in the shape parameter space. Using k-NN, an object is classified by a majority vote of its neighbours, with the object being assigned to the class most common amongst its k nearest neighbours. Third, a SVM constructs a hyperplane or set of hyperplanes in a high-dimensional shape parameter space that can be used for classification, and an optimal separation is achieved by maximising the distance from the hyperplane to the nearest training data point of any class, namely functional margin. In general, we can obtain lower generalisation error of the classifier when there is larger margin.

FS based on statistical hypothesis testing provided guidance of which shape parameters might be used most effectively for classification of the tumour type. In an embodiment, the forward and backward FS can additionally be integrated in the tumour classification procedure. In such an embodiment, in the forward FS method, one shape parameter is selected and then more shape parameters are added iteratively to analyse if the classification performance is improved. If adding a shape parameter can improve the classification performance, the shape parameter is included; otherwise, the shape parameter is excluded. In contrast, the backward FS method starts with using all the shape parameters and gradually removes a shape parameter to test if the classification performance is worsened.

The inventors carried out leave-one-out external cross-validation in classifying GBM and MET by applying three different classification methods, i.e., LDA, k-NN and nonlinear SVM, and two types of shape parameter ranking methods (filter methods and wrapper methods). The statistical assessments were performed to provide the empirical receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC), which is a non-parametric estimate (Metz 1986; Zweig 1993). In addition, the accuracy (Ace), sensitivity (Sens.), specificity (Spec.) of the classifiers with various shape parameters selected were compared.

For the FLDA, the AUC values for the selection of shape parameters when separating GBM from MET varied from 0.86 to 0.91, the highest values being achieved by filter-based FS with three selected shape parameters vCxR, mCvN and stdSI. Compared to the wrapper-based FS, the highest AUC of 0.89 was achieved by using 5 shape parameters. The highest specificity and accuracy were also reached by using three shape parameters selected by the filter-based method, while the lowest specificity and accuracy were found using only one shape parameter. Although only using one shape parameter obtained lowest accuracy and AUC, the values (mean AUC of the Leave-One-Out cross-validation) were over 0.85, which still indicate good performance. Also, we can notice that when we used 10 or 11 features the resubstitution errors were reduced but the cross-validation accuracies were also reduced. Subsequently, we investigated the performance of SVM. The highest AUC of 0.99 was obtained using three shape parameters selected by the filter-based method. Interestingly, both using 5 and 11 shape parameters obtained very high accuracy and AUC (both are 0.99) according to the wrapper-based shape parameter selection.

Compared to the earlier methods of classification between GBM and MET that achieved less than 60% accuracy (Devos et al. 2004; Georgiadis et al. 2008), the present example method obtained a cross- validated accuracy of 86% with 0.91 AUC using the FLDA method and 95% accuracy with 0.99 AUC using the SVM method. A recent study using DTI to discriminate GBM from MET has reported 86% overall classification accuracy (Byrnes et al. 2011). (Blanchet et al. 2011), which discloses a simple model to describe the shape of the tumour in a slice by slice manner and claims a good accuracy of 93.9% has been discussed above; however, the reported accuracy in their paper is the best accuracy of the leave-one-out cross-validation not the mean. Intuitively, a 3D method is complicated and more time-consuming than a 2D method. Compared to Blanchet's 2D method based on stacks of 2D slices, our 3D model is based on a 3D mesh created by the segmentation results, all of the shape parameters being calculated by the voxel numbers of the volume and vertices on the mesh surface, which is computationally more efficient and contains all the information of the segmented tumour region. A further, recent investigation (Chen et al. 2012) classifies GBM and MET using normalised signal intensity, which is an over-simplified model to describe intensity variance within the tumour region. There was no accuracy reported by this study; however, the overall AUC of their ROC analysis is 74.1% and 72.5% for GBM and MET respectively with 0.644 sensitivity and 0.762 specificity.

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 18, 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 66, 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.

References

Awate, Suyash P, Lawrence Win, Paul Yushkevich, Robert T Schultz, and lames C Gee. 2008. "3D Cerebral Cortical Morphometry in Autism: Increased Folding in Children and Adolescents in Frontal, Parietal, and Temporal Lobes." Medical Image Computing and Computer-assisted Intervention: MICCAI ... International Conference on Medical Image Computing and Computer-Assisted Intervention 11 (Pt 1) (lanuary): 559-67. http://www.ncbi.nlm.nih.gov/pubmed/18979791.

Barrick, Thomas R, and Chris A Clark. 2004. "Singularities in Diffusion Tensor Fields and Their Relevance in White Matter Fiber Tractography." Neurolmage 22 (2): 481-491.

doi : 10.1016/j .neuroimage .2004.02.001.

http://www.sciencedirect.eom/science/article/pii/S 1053811904000928.

Blanchet, L, P W T Krooshof, G I Postma, A I Idema, B Goraj, A Heerschap, and L M C Buydens. 2011.

"Discrimination Between Metastasis and Glioblastoma Multiforme Based on Morphometric Analysis of ^'MR Images." American Journal of Neuroradiology 32 (1): 67-73. doi: 10.3174/ajnr.A2269.

http://www.ajnr.Org/content/32/l/67.abstract.

Byrnes, Tiernan J. D., Thomas R. Barrick, B. Anthony Bell, and Chris a. Clark. 2011. "Diffusion Tensor Imaging Discriminates Between Glioblastoma and Cerebral Metastases in Vivo." NMR in Biomedicine 24 (1) (lanuary 23): 54-60. doi: 10.1002/nbm. l555. http://doi.wiley.com/10.1002/nbm.1555.

Chen, X Z, X M Yin, L Ai, Q Chen, S W Li, and I P Dai. 2012. "Differentiation Between Brain

Glioblastoma Multiforme and Solitary Metastasis: Qualitative and Quantitative Analysis Based on Routine MR Imaging." A JNR. American Journal of Neuroradiology 33 (10) (November): 1907-12. doi: 10.3174/ajnr.A3106. http://www.ncbi.nlm.nih.gov/pubmed/22743640. Chernyaev, EV. 1995. "Marching Cubes 33: Construction of Topologically Correct Isosurfaces." Institute for High Energy Physics, Moscow, Russia, Report CN/95-17. 1-8.

http ://citeseerx.ist.psu.edu/viewdoc/download?doi= 10.1.1.56.7139&rep=rep 1 &type=pdf .

Cristianini, N., and J. Shawe-Taylor. 2000. An Introduction to Support Vector Machines and Other Kernel- based Learning Methods . First. Cambridge University Press, http://www.support-vector.net/.

Devos, A, L Lukas, J A K Suykens, L Vanhamme, A R Tate, F A Howe, C Majos, et al. 2004.

"Classification of Brain Tumours Using Short Echo Time 1H MR Spectra." Journal of Magnetic Resonance 170 (1): 164-175. doi: 10.1016/j .jmr.2004.06.010.

http://www.sciencedirect.com/science/article/pii/S 1090780704001818.

Fisher, RA. 1936. "The Use of Multiple Measurements in Taxonomic Problems." Annals of Eugenics 7 (2):

179-188. http://onlinelibraiy.wiley.eom/doi/10. l l l l/j .1469-1809.1936.tb02137.x full.

Frieboes, Hermann B, Xiaoming Zheng, Chung-Ho Sun, Bruce Tromberg, Robert Gatenby, and Vittorio Cristini. 2006. "An Integrated Computational/Experimental Model of Tumor Invasion." Cancer Research 66 (3): 1597-1604. doi: 10.1158/0008-5472.CAN-05-3166.

http://cancerres.aacrjournals.Org/content/66/3/1597.abstract.

Friedman, JH, JL Bentley, and RA Finkel. 1977. "An Algorithm for Finding Best Matches in Logarithmic Expected Time." A CM Transactions on Mathematical Software (TOMS) 1549 (July): 209-226.

http : //dl .acm .org/citation . cfm ?id=355745.

Georgiadis, Pantelis, Dionisis Cavouras, Ioannis Kalatzis, Antonis Daskalakis, George C Kagadis, Koralia Sifaki, Menelaos Malamas, George Nikiforidis, and Ekaterini Solomou. 2008. "Improving Brain Tumor Characterization on MRI by Probabilistic Neural Networks and Non-linear Transformation of Textural Features." Computer Methods and Programs in Biomedicine 89 (1): 24-32.

doi: 10.1016/j .cmpb.2007.10.007.

http://www.sciencedirect.com/science/article/pii/S0169260707002404.

Gong, Faming, and Xinglong Zhao. 2010. "Three-Dimensional Reconstruction of Medical Image Based on Improved Marching Cubes Algorithm." 2010 International Conference on Machine Vision and Human- machine Interface (1987): 608-611. doi: 10.1109/MVHI.2010.45.

http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5532558. Ter Haar Romeny, Bart M. 2011. "Biomedical Image Processing." In , ed. Thomas M Desemo, 183-187. Springer.

Haury, Anne-Claire, Pierre Gestraud, and Jean-Philippe Vert. 2011. "The Influence of Feature Selection Methods on Accuracy, Stability and Interpretability of Molecular Signatures." PloS One 6 (12) (January): e28210. doi: 10.1371/journal.pone.0028210.

http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3244389&tool=pmcentrez&rendertype=abs tract.

Koenderink, Jan J, and Andrea J van Doom. 1992. "Surface Shape and Curvature Scales." Image and Vision Computing 10 (8): 557-564. doi: 10.1016/0262-8856(92)90076-F.

http://www.sciencedirect.com/science/article/pii/026288569290076F.

Lachenbruch, Peter A., and M. Goldstein. 1979. "Discriminant Analysis. " Biometrics: 69-85.

Law, Meng, Soonmee Cha, Edmond A Knopp, Glyn Johnson, John Arnett, and Andrew W Litt. 2002.

"High-Grade Gliomas and Solitary Metastases: Differentiation by Using Perfusion and Proton

Spectroscopic MR Imaging l ." i?ad/o/ogy 222 (3): 715-721. doi: 10.1148/radiol.2223010558.

http : //radiology . rsna. org/content/222/3/715. abstract.

Levner, Ilya. 2005. "Feature Selection and Nearest Centroid Classification for Protein Mass Spectrometry." BMC Bioinformatics 6 (January): 68. doi: 10.1186/1471-2105-6-68.

http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1274262&tool=pmcentrez&rendertype=abs tract.

Lewiner, Thomas, and Hello Lopes. 2003. "Efficient Implementation of Marching Cubes' Cases with

Topological Guarantees." Journal of Graphics Tools (January 2013): 37-41.

http://www.tandfonline.com/doi/full/10.1080/10867651.2003.10487582.

Lin, Qingfen. 2003. "Enhancement, {D}etection, and {Visualization of 3 {D} {V}olume {D}ata". SE-581 83 Linkoping, Sweden: Dept. EE.

Liu, Huan, and Lei Yu. 2005. "Toward Integrating Feature Selection Algorithms for Classification and

Clustering." Knowledge and Data Engineering, IEEE Transactions On 17 (4): 491-502.

http : //ieeexplore . ieee .org/xpls/abs_all j sp ?arnumber= 1401889. Lorensen, WE, and HE Cline. 1987. "Marching Cubes: A High Resolution 3D Surface Construction Algorithm." A CM Siggraph Computer Graphics 21 (4): 163-169.

htt : //dl .acm .org/citation . cfm ?id=37422.

Metz, CE. 1986. "ROC Methodology in Radiologic Imaging." Investigative Radiology 21 : 720-733.

http://joumals.lww.com/investigativeradiology/Abstract/1986/09000/ROC_Memodology_in_Radiologi c_Imaging .9.aspx.

Newman, Timothy S., and Hong Yi. 2006. "A Survey of the Marching Cubes Algorithm." Computers & Graphics 30 (5) (October): 854-879. doi: 10.1016/j .cag.2006.07.021.

http://linkinghub.elsevier.com/retrieve/pii/S0097849306001336.

Nielson, GM. 2003. "On Marching Cubes." Visualization and Computer Graphics, IEEE Transactions On 9 (3): 283-297. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1207437.

Opstad, Kirstie S, Mary M Murphy, Peter R Wilkins, B Anthony Bell, John R Griffiths, and Franklyn A Howe. 2004. "Differentiation of Metastases from High-grade Gliomas Using Short Echo Time 1H Spectroscopy." Journal of Magnetic Resonance Imaging 20 (2): 187-192. doi: 10.1002/jmri.20093.

Saeys, Yvan, Inaki Inza, and Pedro Larranaga. 2007. "A Review of Feature Selection Techniques in

Bioinformatics." Bioinformatics (Oxford, England) 23 (19) (October 1): 2507-17.

doi: 10.1093/bioinformatics/btm344. http://www.ncbi.nlm.nih.gov/pubmed/17720704.

Senger, Donna, J Gregory Cairncross, and Peter A J Forsyth. 2003. "Long-Term Survivors of Glioblastoma:

Statistical Aberration or Important Unrecognized Molecular Subtype?" The Cancer Journal 9 (3): 214- 221.

http://journals.lww. com/journalppo/Fulltext/2003/05000/Long_Term_Survivors_of_Glioblastoma Sta tistical.9.aspx.

Soffietti, Riccardo, Roberta Ruda, and Roberto Mutani. 2002. "Management of Brain Metastases." Journal of Neurology 249 (10): 1357-1369. doi: 10.1007/s00415-002-0870-6. http://dx.doi.org/10.1007/s00415- 002-0870-6.

Suykens, JAK, and J Vandewalle. 1999. "Least Squares Support Vector Machine Classifiers." Neural

Processing Letters (9): 293-300. http://link.springer.com/article/10.1023/A%3A1018628609742. Tsuchiya, K, A Fujikawa, M Nakajima, and K Honya. 2005. "Differentiation Between Solitary Brain Metastasis and High-grade Glioma by Diffusion Tensor Imaging." British Journal of Radiology 78 (930): 533-537. doi: 10.1259/bjr/68749637. http://bjr.birjournals.org/content/78/930/533.abstract.

Vellido, A, E Romero, M Julia-Sape, C Majos, A Moreno-Torres, J Pujol, and C Anis. 2012. "Robust Discrimination of Glioblastomas from Metastatic Brain Tumors on the Basis of Single-voxel (1)H MRS." NMR in Biomedicine 25 (6) (June): 819-28. doi: 10.1002/nbm. l797.

http : //www .ncbi .nlm .nih .gov/pubmed/22081447.

Zweig, H. 1993. "Receiver-Operating Clinical Medicine ( ROC ) Plots: A Fundamental Evaluation Tool In" 39 (4): 561-577.

Claims

1. 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;

identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels;

processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape;

determining the tumour type based on the one or more calculated shape parameters.

2. A method according to claim 1, further comprising:

identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein:

one or more of the shape parameters is/are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches.

3. A method according to claim 2, wherein the one or more of the shape parameters that are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches comprise(s) one, more than one in any combination, or all of the following:

average Shape Index, avgSI;

standard deviation of the Shape Index, stdSI;

mean Shape Index, mSI;

standard deviation of the Curvedness, stdCv;

mean Curvedness, mCvN; and

the mean curvature L2-norm, mCN, which is defined as follows:

m N = _ieS H² (i ds ⁵ .

wherein Ή is the mean curvature of each surface patch and S is the tumour surface.

4. A method according to claim 3, wherein the one or more of the shape parameters that are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches comprise(s) one, more than one in any combination, or all of the following:

average Shape Index, avgSI;

standard deviation of the Shape Index, stdSI;

mean Curvedness, mCvN; and the mean curvature L2-norm, mCN, which is defined as follows:

m N = _ieS H² (i ds ⁵ .

wherein Ή is the mean curvature of each surface patch and S is the tumour surface; and

the determining of the tumour type is further based on either or both of the volume ratio of the tumour with respect to the minimum convex hull, vCxR, and the volume ratio of the tumour with respect to the minimum cuboid bounding box, vCR.

5. A method according to claim 4, wherein the one or more of the shape parameters that are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches comprises:

avgSI or stdSI but not both.

6. A method according to claim 4 or 5, wherein the one or more of the shape parameters that are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches comprises:

mCN or mCvN but not both.

7. A method according to any of claims 2-6, wherein a combination of two or more local three- dimensional shape characteristics is determined for each of the plurality of surface patches, and at least one of the shape parameters is derived based on a distribution of the combination of two or more local three- dimensional shape characteristics over the plurality of surface patches.

8. A method according to claim 7, wherein the combination of two or more local three-dimensional shape characteristics comprises the Shape Index, SI, and the Curvedness, Cv.

9. A method according to claim 7 or 8, wherein at least one of the shape parameters is derived based on a density distribution or histogram representing the distribution of the combination of two or more local three-dimensional shape characteristics over the plurality of surface patches.

10. A method according to any of the preceding claims, wherein the one or more shape parameters comprises one, more than one in any combination, or all of the following:

the standard deviation of the Normalized Radial Length, stdCD, the radial length being defined as the distance from the centroid of the volume inside the tumour surface and the normalised radial length being defined as the distance between the centroid and the tumour surface normalized by the maximum value of the radial length; the volume ratio of the tumour with respect to the minimum cuboid bounding box, vCR; the volume ratio of the tumour with respect to the minimum convex hull, vCxR.

11. A method according to claim 10, wherein the one or more shape parameters comprises:

vCxR or vCR but not both.

12. A method according to any of the preceding claims, further comprising:

identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein:

the one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape comprise mCvN, stdSI and vCxR.

13. A method according to any of the preceding claims, wherein the determination of the tumour type comprises distinguishing between a glioblastoma multiforme, GBM, and a solitary metastasis, MET.

14. A method according to any of the preceding claims, wherein the MRI data comprises Diffusion Tensor Imaging, DTI, data and the segmenting of the MRI data to identify tumour voxels and non-tumour voxels is performed based on decomposing the diffusion tensor into isotropic, p, and anisotropic, q, components and comparing the values of p and q for each voxel with predetermined maps defining regions in p, q space that correspond to different tissue types.

15. A method according to any of the preceding claims, wherein the one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape are each representative of a three-dimensional characteristic of the whole surface shape.

16. 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;

identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels;

processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape; determining the tumour type based on the one or more calculated shape parameters.

17. A data analysis unit according to claim 16, wherein the data processing unit is further configured to perform the following steps:

identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein:

one or more of the shape parameters is/are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches.

18. 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;

identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels; and

processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape;

determining the tumour type based on the one or more calculated shape parameters.

19. An MRI machine according to claim 18, wherein the data processing unit is further configured to perform the following steps:

identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein:

one or more of the shape parameters is/are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches.

20. 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;

identifying a tumour surface representing a boundary between the tumour voxels and the non-tumour voxels;

processing the identified tumour surface to calculate one or more shape parameters that are each representative of a three-dimensional characteristic of the surface shape;

determining the tumour type based on the one or more calculated shape parameters.

21. A computer program according to claim 20, further comprising:

identifying a plurality of surface patches representing different portions of the tumour surface; and determining one or more local three-dimensional shape characteristics of each of the plurality of surface patches, wherein:

one or more of the shape parameters is/are derived based on a distribution of the one or more local shape characteristics over the plurality of surface patches.

22. A method of analysing MRI data of a tumour to determine the tumour type substantially as hereinbefore described with reference to and/or as illustrated in the accompanying figures.

23. A data analysis unit or MRI machine configured and arranged to operate substantially as hereinbefore described with reference to and/or as illustrated in the accompanying figures.