DESCRIPTION

TECHNICAL FIELD

The present invention relates to a computerized method for classifying populations of objects described by n-dimensional matrices. In details, the invention relates to a method for classifying populations of objects described by n-dimensional matrices of scalar values measured within a space wherein the variability range of the scalar values is comparable with the noise of the measurement and the function representing the scalar values within said space, or at least in separable subsets, is smooth.

The invention is preferably and advantageously applied in processes analyzing 3D images, characterized by a low contrast and/ or a low resolution, such as PET (Positron Emission Tomography) images or SPECT (Single Photon Emission Tomography) images.

In particular the invention is preferably and advantageously applied as a part of a process analyzing and classifying PET or SPECT images of brains, intended to help the physicians in the diagnosis of neurodegenerative diseases such as Alzheimer's disease, Parkinson disease or Dementia with Lewy bodies.

PRIOR ART

Mathematical methods for processing scalar 3D matrices are widely used in several fields, for example for processing 3D images, which simply are scalar 3D matrices whose elements are values of pixels of the image.

In the medical field, in particular, PET and SPECT produce 3D images of the observed tissues; the processing of such images allows the physicians to have information necessary to make the diagnosis.

Alzheimer's disease is a disease difficult to be diagnosed without the help of invasive techniques and the use of several biomarkers. In particular the certification of a death caused by said disease can be confirmed only with a histopathologic examination of the brain tissue, namely cerebral tissue, taken post mortem, on which the analysis

l shows the presence of neuritic plaques, composed of depositions of amyloid-β peptides (particularly Αβ 1-40 and Αβ 1-42).

Parkinson disease and the Dementia with Lewy bodies on the contrary are characterized by a dopamine deficiency, caused by the death of neurons in the substantia nigra. The reason of such degeneration of the neurons is unknown, but it seems to be associated to the presence of aggregates of the a-synuclein protein.

The presence in the brain of proteins that are precursors of such diseases currently is considered as a considerable risk factor and therefore investigation techniques have been developed for assessing the pathologic protein load. By injecting in a patient radiopharmaceuticals able to be bound to amyloid-β or a-synuclein proteins , it is possible to obtain PET or SPECT images that map the presence of such proteins in the brain and therefore they allow the physician to make a diagnosis. In order to give to the physician useful information, PET apparatuses are known able to process PET images for assessing amyloid-β protein concentrations. Software used by such apparatuses are proprietary software, depending on the particular apparatus, that require very high computational resources and that suffer from instability, since they are based on the quantification and the ratio of PET counts in relatively limited cerebral regions, whose identification (segmentation) is particularly difficult.

OBJECTS AND SUMMARY OF THE INVENTION

It is the object of the present invention to overcome the prior art drawbacks.

In particular, it is the object of the present invention to provide a low-cost computational method that allows populations described by a 3D matrix to be classified without the need of a dedicated computer.

In particular it is an object of the present invention to provide a low-cost computational method that provides to a physician an index useful for the diagnosis of neurodegenerative diseases, in particular a method allowing distinguishing between subjects at risk or not at risk of neurodegenerative diseases such as the Alzheimer's disease or Parkinson disease.

These and other objects of the present invention are achieved by a method embodying the characteristics of claim 1. Such method provides to compute, for each image - or in general n-dimensional matrix - two different metrics. The image, or matrix, therefore is classified as a function of its distance from a joined metric, determined on the basis of the first two metrics computed for a statistical sample of images or matrices.

Unlike known methods classifying PET images, that require a considerable computational power for segmenting the image into very small portions, the method analyzes geometrical or scalar macro-features of the analyzed image/matrix; therefore the method requires a smaller computational power and it can be implemented on common devices and not on proprietary devices.

Moreover the suggested method has interesting scalability characteristics, since the computation of the two metrics can be carried out in parallel by several distinct machines or central processing units.

Further advantageous characteristics of the present invention will be more clear from the following description and from the annexed claims, which are an integral part of the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will be described below with reference to non limitative examples, provided by way of example and not as a limitation in the annexed drawings. These drawings show different aspects and embodiments of the present invention and, where appropriate, reference numerals showing like structures, components, materials and/ or elements in different figures are denoted by like reference numerals.

Figure 1 is a flow chart of a method according to the invention;

Figures 2a and 2b are each one two PET images, before (top) and after (bottom) a normalization process;

Figure 3 is a flow chart of the process for computing a first metric of the method of figure 1;

Figure 4 is some intermediate results of the computation of the first metric for a PET image;

Figure 5 is a flow chart of the process for computing a second metric of the method of figure 1;

Figure 6 is a flow chart for the computation of a joined metric to be used in the classification of a PET image;

Figure 7 is the experimental results in the classification of Amyloid-PET images. DETAILED DESCRIPTION OF THE INVENTION

While the invention is susceptible of various modifications and alternative forms, some non limitative embodiments, provided by way of example, are described here below in details.

It should be understood, however, that there is no intention to limit the invention to the specific embodiments disclosed, but, on the contrary, the intention of the invention is to cover all modifications, alternative forms and equivalents falling within the scope of the invention as defined in the claims.

Therefore in the description below the use of "for example", "etc", "or" indicates non-exclusive alternatives without limitation unless otherwise defined; the use of "also" means "among which, but not limited to" unless otherwise defined; the use of "includes/comprises" means "includes/ comprises, but not limited to," unless otherwise defined.

The term amyloid-PET refers to a PET obtained by injecting in the patient radiopharmaceuticals or molecules able to be bound to amyloid proteins.

The term voxel refers to an element of elementary volume comprising a three- dimensional volume wherein detections are made. By way of example, a cerebral PET image typically is composed of 256 x 256 x 128 voxels.

In the description below, the method according to the present invention is disclosed with reference to an application in the medical field, such application however has not to be intended as a limitation of the invention, which on the contrary relates to an analysis method that is generally applied to sets of homogeneous data that can be represented as n-dimensional matrices containing a scalar datum. The method extracts geometric features based on the definition of suitable hypersurfaces and volumes, and it uses such information in order to associate to each matrix of the set a metric derived from the association of at least two independent metrics that can be used as a discriminant measurement in the classification of the matrix.

In the example below, described with reference to the flow chart of figure 1, the method is directed to classify brain images - for example PET or SPECT images - into groups of images related to potentially healthy subjects and images related to potentially pathological subjects. Such classification is not a diagnosis, but it can be useful both for merely image reorganization purposes, and as an index for a physician that then has to make the diagnosis.

The method starts at step 100 of figure 1, wherein a PET image of a brain is acquired (for example of the type shown in figure 2a). Figure 2 shows two different subjects (2a and 2b), acquired according to two PET imaging protocols and therefore with a different signal-to-noise ratio. The top portion of the figures (denoted by the reference 200) is a coronal section of the image directly acquired by the scanner, while the lower part of the figures (denoted by the reference 201) shows a coronal section of the result of a first standardization process. Figures 2a and 2b therefore allow the strenght of the method to be appreciated also in presence of images with a non similar quality.

In one embodiment, the acquisition step comprises the process of creating the image by a PET machine, as an alternative such step is limited to the acquisition of an electronic file generated by a PET machine that generates an output image in the form of an electronic file (figures 2a and 2b, top portion). As it is known, PET images are 3-dimensional matrices acquired by injecting radiopharmaceuticals marked with radioactive elements that bind in preferential sites useful for determing brain pathologies. The matrices contain scalar values, that represent the counts due to the radiopharmaceutical uptake in the individual volume element (voxel). In the example described below, PET images make reference to a patient injected with a radiopharmaceutical bound to the amyloid protein. Such images therefore are useful for the early diagnosis of the Alzheimer's disease.

Once the image is acquired and pre-processed to allow data coming from different scanners and acquisition protocols to be compared, the method of figure 1 provides (step 101) to calculate a first metric based on geometric features of the image, to calculate a second metric (step 102) based only on intensity values of the image, and to calculate then a joined metric (step 103) that allows the image to be discriminated/ classified as a function of the calculated values of the first and second metrics.

Now with reference to the details of the method, as said above the first metric considers the geometric features of the image, therefore in one embodiment, the process computing the first metric (shown with reference to figure 3) provides (step 301) to carry out a spatial normalization of the image such to bring the acquired image within predetermined spatial features. Such step is optional and it is not necessary if data have been already acquired under isometric conditions.

In one embodiment, the spatial normalization step for a PET image of a brain provides the following steps:

a) Correcting anomalous intensities: voxels having anomalous values are corrected, for example voxels higher than a given intensity are brought to a maximum accepted intensity value

b) Filtering the image for adapting the image to a standard model for example in the case of brain images, the standard model may be a standard MRI template (CBM152, Montreal Neurological Institute [1]). The registration is multimodal (PET -> MRI) and it uses specific filters for making the reference MRI compatible with the PET. Since information contained in resonance images and PET images are physically different, it has been necessary to develop a multimodal normalization and registration procedure in order to align the corresponding anatomical structures in the different considered subjects. Such procedure uses both known and original methods. Particularly the use of decomposition into "curvelet" functions for the standard MRI template is mentioned, without this the registration between MRI and PET would not be strong [4] [5].

Figures 2b and 2a in the lower part show the result after the spatial normalization process. After the spatial normalization, the method provides to define (step 302) a region of interest (ROI) of the image. Such region of interest can coincide with the whole image or can comprise only one part thereof.

In the application example of the analysis of PET images for the diagnosis of the Alzheimer's desease, this step consists in extracting the cortical surface and the surface of the ventricles of the brain. In order to accurately define these surfaces, ROIs pre-segmented on the MRI template are preferably used, whose identification has been already clinically validated [6]. These pre-segmented ROIs however are not adapted to the different morphologies of the subjects, therefore a method for adapting them is necessary. The adaptation preferably occurs by using the distance transform, that is a function derived from the binary image of the ROI. Such function generates a 3D map that is modulated (namely convoluted) with the PET, such to adapt the contour (boundary) of the pre-segmented ROI on the counts of the PET. A similar treatment is also suggested on ventricle volumes, such that the cortical surface and the contour of the ventricle volumes represent the limit surfaces of the ROI. Morvoer, depending on the application and on the data, such step can be further improved by defining a region of interest (ROI), that in the amyloid-PET application is exemplified by the selection only of the lobes and sub-cortical structures, excluding the cerebellum and the brain stem.

Once defined the region of interest, the method provides to normalize (step 303) the scalar values of the image. In the example shown herein, such step is carried out by bringing the intensity values of the voxels of the image back into the range [0,1] by evaluating the quantiles and the residual counts in the ventricle volumes and within said limit surfaces. In practice, such step considers (even if in a perfunctory manner) the problem of the partial volume in PET counts. In the cerebrospinal fluid regions uptake should not occur, but in the images there are counts also in the cerebrospinal fluid due to the low spatial resolution and to the low signal-to-noise ratio related to the imaging technique (partial volume effect). In the normalization of intensities the level of the counts within the ventricles (determined by previous ROIs) and those delimited within the cortical surface will be considered. The calulation of the quantiles on the distribution of the intensities delimited in this manner will set to 0 the level of intensities found in the ventricles, and to 1 the 99% quantile of the distribution. The result is a change in the scale of intensities with the exclusion (namely intensities set at 0) of some regions (ventricles) and the maximum set at 99% of the total distribution of the counts.

This normalization step is carried out only once on the large ROI described above. If a smaller sub-region is analyzed (e.g. only the front right lobe), the normalization of the intensities would remain the same. The reason clearly is to scale the total counts such to make data coming from different objects and scanners comparable. For this reason it is preferred to have a relatively extended volume and comprising cerebrospinal fluid regions.

Once the normalization of scalar values is carried out, the method provides to calculcate (step 304) N isosurfaces of intensity Ki in (0,1) in the ROI, that is surfaces of the cortical region of the brain that have the same radiation intensity and therefore the same voxel value.

For each isosurface Ki, the relevant area ai and the volume Vi defined as the set of voxels x whose k_x is k_x<= Ki, namely the volume whose boundary is the isosurface is computed (step 305).

Then (step 306) some geometric features of the surfaces ai and of volumes Vi, are computed such to define a first metric for the classification of the 3D matrix.

In a particularly advantageous embodiment for the classification of PET images, the evaluated features are the equivalent radii ¾ and rvi for surfaces ai and for the volumes Vi respectively. Equivalent radii ¾ and iVi mean the radii of spheres having surface equal to ai and volume equal to Vi respectively.

Thus, for each isosurface Ki a pair of values ¾ and r^ are obtained.

By arranging such values in the Cartesian plane having r^ in abscissa and r _a i in ordinate, it is possible to obtain a diagram such as that of figure 4 and to define (step 308) a function that expresses the geometric features of the isosurfaces Ki of the considered 3D matrix. The function F=r _a i(rvi) for example is obtained by interpolation (linear, polynomial, spline or other type) of the points (r ai,f i)^ as an alternative the function rai(rvi) can be obtained as a series of step functions with different values in different ranges of equivalent radii iVi.

In the preferred embodiment a large number of isosurfaces is selected, particularly N higher than or equal to 16, the function F=r _a i(rvi) therefore is obtained by linear interpolation.

The integral of the function defined anyhow, provides (step 309) a first metric ml for the classification of the 3D matrix.

In a preferred embodiment before defining the function F=r _a i(rvi), the values ¾ are normalized (step 307) with respect to the first and last point, that is the points with minimum and maximum rvi (limit surfaces). This step allows the isosurface curves to be made comparable, that is the curves F calculcated for the different isosurfaces Ki. The limit surfaces for all the subjects are the cortical surface and the ventricle surface. By establishing that the point (r _a i,rvi) i=N is the same for all the subjects is the same as considering the morphology of the cortical limit surface therebetween to be uniform. Likewise we can establish the equivalent condition with the minimum limit surface, namely the ventricle one. Now the difference between the subjects only depends on the trend of the function F and not on the starting points thereof. The integral of F normalized in this manner is the first metric ml*.

As mentioned above, the method provides to compute a second metric. Such computation can be advantageously carried out in parallel to the computation of the first metric to reduce the computation time, however it can be carried out before or after computing the first metric.

In the example below, for computing the second metric, only scalar information of the image are used, preferably by supervised partitioning methods (clustering).

In the example described herein with reference to figure 5, the intensity components of the image in the volume of the ROI (step 601) are selected. Such intensity components, that actually are scalar values of the 3D matrix composing the image, then are partitioned into classes (step 602) by a partitioning method, preferably k- means method with two classes, for example as described in the documents [2] or [3] mentioned at the end of the description. k-means method with two classes, for example, allows two classes CI and C2 to be identified with nl and n2 voxel respectively, and whose mean intensities (that is the average of the intensities of the voxels belonging to the class) are Intl and Int2. The method therefore provides to compute (step 603) the second metric as the ratio between the intensities of the two centroids: m2= Intl/Int2.

In a preferred embodiment, the second metric is computed by normalizing the ratio between the classes for the relevant median intensities. In this embodiment, therefore a second metric m2*= (nl*Intl) / (n2*Int2) is provided.

Finally the method provides to determine a joined metric that, on the basis of the two metrics calculated above, allows the image to be classified.

The joined metric is empirically calculated as described below with reference to figure 6, by taking a statistical sample of images (step 701) whose classification is known, and by computing the metrics ml and m2 or ml* and m2* for each sample image (step 702). Each sample image is thus identified by a point, e.g. (ml*, m2*) of the Cartesian plane having the two metrics as abscissa and ordinate.

The joined metric therefore is identified (step 703) as the function that allows the sample images in the Cartesian plane to be divided into homogenous groups.

For example figure 7 shows the result of the process computing the joined metric for the example of PET images described above. The joined metric, in this case, has been computed by taking Amyloid-PET images of 150 subjects coming from the database of the Alzheimer's Disease Neuroimaging Initiative (ADNI). After having computed the metrics ml* and m2* for each image, the points (ml* and m2*) have been arranged on the Cartesian plane. The anaylis of the results shows that the images related to subjects considered as clinically not ill (that is negative, term "neg" in the legend) are positioned in the Cartesian plane in a manner well distinct from ill subjects (that is positive ones, term "pos" in the legend).

Figure 7 further distinguishes the subjectes classified in open condition (namely the subjects whose clinical diagnosis was known to the proposers before the analysis) and subjects classified in blind condition, namely the subjects whose clinical diagnosis (positive or negative with respect to the amyloid load) has been known to the proposers only after the analysis and after reading the images by two independent skilled physicians. Therefore the method allows a cut-off line to be determined by means of which the subjects are classified. Such cut-off line is a joined metric, that is an index, that can be used for the future classification of Amyloid-PET images.

From the numerical point of view, it is appreciated that the method applied to the classification of Amyloid-PET images, gives a reproducibility, with respect to the visual assessment of PET images by skilled physicians, equal to 95% [92-99 with a confidence level (CL) of 95%] (CL is the discriminating power of the negative/ positive subjects in the examination measured with the area under the ROC curve (receiver operating characteristic)).

In the light of the description above it is clear how the method of classification of n- dimensional matrices allows the above objects to be achieved, allowing a very accurate classification, as it is the case of Amyloid-PET images, with a reduced computational effort, and however, without the need of dedicated hardware.

It is also clear that the embodiment described above has not to be intended as a limitation of the present invention and many variants can be made by the person skilled in the art without for this reason departing from the scope of protection as it results from the annexed claims.

In particular, it is clear that the method mentioned above can be carried out as a program for a computer able to execute the method steps once it is run on a computer, both a personal computer, a laptop, a mobile device, a smart phone or another electronic device equipped with an electronic processing unit able to process data.

It is also clear that the method described above with reference to PET image classification can be equally reproduced for classifying images of other type, or more in general, populations collected in n-dimensional matrices. Only by way of example, the n-dimensional matrix instead of being a PET image may be a matrix containing values of electromagnetic fields. The method, in this case, would allow areas of the matrix space to be classified into areas with high electric field and other areas with a limited electric field.

Also in this case, for the classification of the n-dimensional matrix it would be possible to compute a first and a second metric and to define a joined metric by means of which the matrix can be classified depending on the specific metrics calculated.

It is therefore clear that, besides the specific example described above, the invention relates to a computerized method for classifying populations of objects described by n-dimensional matrices. The method comprises the steps of:

a) Electronically acquiring a n-dimensional matrix containing homogeneous scalar values;

b) Computing a first metric (e.g. ml, ml*), said computation comprising the steps of:

i. Defining a region of interest as a subset of the scalar values of the n- dimensional matrix;

ii. Finding N isosurfaces of intensity in the region of interest;

iii. For each isosurface of intensity, calculating the area (ai) of the isosurface and of the volume (vi) underneath said isosurface;

iv. For each isosurface of intensity calculating a point (r _a i , rvi) in the Cartesian plane Tai , Γνί , wherein

rai is the equivalent radius of a sphere having the surface equal to the area of the isosurface, and

rvi is the equivalent radius of a sphere having the volume equal to the volume underneath the isosurface;

v. Defining a function rai(rvi) that passes by or that is obtained by interpolation of the N points (r _a i , rvi) calculated for each isosurface. vi. Defining the first metric as the integral of the function r _a i(rvi). c) Computing a second metric (m2, m2*), said computation of the second metric (m2, m2*) comprising the steps of:

i. Selecting a set of scalar values of the region of interest; ii. Partitioning the set of scalar values in a number m of classes, m being an integer greater than 1,

iii. Defining the centroids of the m classes in which the set of scalar values of the matrix has been partitioned,

iv. Defining the second metric as the m-1 dimensional vector whose elements are ratios between scalar values of the centroids of two of said m classes multiplied by the number of elements belonging to the two classes.

d) Defining a joined metric as a function able to separate, in the two-dimensional space of the first and second metrics, the matrices of the sample of a first class from the other matrices of the sample.

e) Classifying the n-dimensional matrix acquired at step a) depending on its distance from the joined metric.

Advantageously before defining the function rai(rvi), the method provides to proceed to normalize the points (r _a i, rvi) such that the point (¾, rvi) with maximum rvi , coincides with a reference point (r _ama x, rvmax) of a template. More preferably the normalization provides to cause also the point (¾, rvi) with minimum r^ to coincide with a reference point ( in, rvmin) of a template. This allows the scalar matrix to be classified (for example the image of a brain) by comparig it with homogeneous matrices (e.g. images of brains rescaled on the same template by means of rotations and/or translations and/or trapezoidal or "shear" deformations and/or scale changes).

Advantageously, then the method can provide to compute the second metric by normalizing the ratio between the number of scalar values belonging to the classes determined for the relevant median intensities.

Then preferably the determination of the joined metric comprises the steps of:

i. Acquiring a sample of matrices already classified and homogeneous to the n-dimensional matrix to be classified;

ii. For each matrix of the sample computing the first metric and the second metric according to the steps mentioned above as b) and c) of the method, such to associate a point in the two-dimensional space of the first and second metrics to each matrix of the sample; iii. Defining the joined metric as a function able to separate, in the two- dimensional space of the first and second metrics, the matrices of the sample of a class from the matrices of the sample of other classes.

A characteristic of the method is to capture small differences in the global distribution of information into the image/ matrix.

Therefore the method is particularly advantageous if applied for classifying n- dimensional matrices characterized by continuous distribution and by a variability range reduced and comparable with the noise.

Therefore the method can be advantageously applied if the n-dimensional matrix to be classified corresponds to an image of a brain, particularly a PET image, where the scalar values of the matrix represent intensity values of emissions produced by radiopharmaceuticals.

In the case of the classification of a brain image, after having electronically acquired the image and before calculating the first and second metrics, it is possible to advantageously proceed with a spatial normalization of the image, by filtering the acquired image such to adapt it to a standard brain model, particularly the ICBM152 model, by Montreal Neurological Institute.

Still in the case of the classification of a brain image, in the computation of the first metric the region of interest can be advantageously selected as the region comprised between the cortical surface and the contour of the ventricle volumes of the brain, and particularly it is composed of the lobes and of the sub-cortical structures excluding the cerebellum and brain stem. Obviously other regions may be selected they being advantageous depending on the desired classification.

In order to give to the physician information for making a diagnosis, it is possible to display the position of the acquired PET image in the plane ml,m2 or ml*,m2* together with the other sample images previously classified and/ or with the joined metric. As an alternative it is possible to provide an output datum from the computer, for example in the numerical form, indicative of the distance (in the plane ml,m2 or ml*,m2*) of the acquired image from the joined metric.

References

[1] Mazziotta, J. C, Toga, A. W., Evans, A., Fox, P., & Lancaster, J.(1995). A probabilistic atlas of the human brain: theory and rationalefor its development. The International Consortium for Brain Mapping(ICBM). Neuro Image, 2(2), 89-101.

[2] Seber, G. A. F. (2008) Multivariate Distributions, in MultivariateObservations, John Wiley & Sons, Inc., Hoboken, NJ, USA.

[3] Mucha, H.-J. (1986), Spath, H.: Cluster dissection and analysis: theory, FORTRAN programs,examples. (Translator: Johannes Goldschmidt.) EllisHorwood Ltd Wiley, Chichester 1985. 226 pp. £25. Biom. J., 28: 182. doi:10.1002/bimj.4710280207

[4] M. N. Safran; M. Freiman; M. Werman; L. Joskowicz: Curvelet-based sampling for accurate and efficient multimodal image registration. Proc. SPIE 7259, Medical Imaging 2009: Image Processing, 72590M (27 March 2009); doi: 10.1117/12.810695

[5] Paolo Bosco, Andrea Chincarini, Mario Esposito, Luca Rei, Cesare Ovando: MRI and PET combined analysis of the medial temporal lobe. Medical Measurements and Applications Proceedings (MeMeA), 2011 IEEE International Workshop 01/2011; DOI: 10.1109/MeMeA.2011.5966721

[6] VS Fonov, AC Evans, K Botteron, CR Almli, RC McKinstry, DL Collins and BDCG. Unbiased average age-appropriate atlases for pediatric studies. NeuroImage,Volume 54, Issue 1, January 2011, ISSN 1053-8119, DOI: 10.1016/j.neuroimage.2010.07.033