Non-invasive joint evaluation

[0001] This application claims the benefit of U.S. Provisional Application No. 60/584,663, filed Jul. 1, 2004, which is herein incorporated by reference.

[0002] This invention relates to methods and apparatus for evaluating the condition of one or more joints of a patient using a non-invasive data acquisition technique, such as Magnetic Resonance Imaging (MRI).

[0003] Osteoarthritis (OA) is a prevalent disease characterized mainly by cartilage degradation that is clinically reflected by a gradual development of joint pain, stiffness, and loss of motion. Osteoarthritis is extremely frequent in the general population, and it is estimated that its radiological prevalence is close to 50% overall. This figure is even higher in the elderly, with as much as 75% of the population between ages of 55 and 64 exhibiting some degree of radiological osteoarthritis in one or more joints. Although this disease is often benign, severe degenerative changes may cause serious disability.

[0004] The early stages of osteoarthritis are relatively asymptomatic. Patients then begin to experience mild symptoms, such as morning stiffness, pain with activity, and mild swelling. In later stages, the frequency and intensity of pain tends to increase, and there can be a progressive loss of joint function.

[0005] Current non-surgical treatment regimens focus on maintaining function and decreasing symptoms. These can include the use of anti-inflammatory medications, exercise programs, weight control, and occasionally steroid injections. Beginning these regimens at earlier stages may delay the onset of more serious symptoms, and can in some cases even postpone or eliminate the need for surgery.

[0006] Promising research is also under way to find a therapeutic agent that will slow or stop the progression of the disease. But it is unlikely that this research will yield an agent that will actually reverse the progression of the disease. Even if current efforts are successful in identifying suitable therapeutic agents, therefore, the overall comfort levels and joint function levels achievable with these agents will most likely still depend on early and accurate diagnosis.

[0007] In one general aspect, the invention features a musculoskeletal imaging system that includes a source of feature data extracted from imaging data resulting from imaging acquisitions from joints of different individuals affected by different diseases, and this feature data includes disease characteristic categorization information for a plurality of disease categories. A comparison module is operative to compare patient imaging data resulting from an imaging acquisition from a joint of a patient with the feature data. The comparison module is also operative to provide at least one categorization indicator for the patient imaging data that indicates a correspondence between spatial information in the patient imaging data and the disease categories for which there is extracted categorization information in the feature data.

[0008] In preferred embodiments, The source of feature data can be operative to provide feature data for at least some individuals having healthy cartilage. The source of feature data can be operative to provide feature data that expresses information for a plurality of patients that has been normalized and aggregated. The source of feature data can be operative to provide feature data that expresses information for a plurality of patients that has been normalized, aggregated, and compressed. The source of feature data can be operative to provide feature data that expresses information for a plurality of patients that has been normalized and averaged. The source of feature data can be operative to provide feature data that expresses information for a plurality of segmented patient scans that have been normalized and aggregated. The source of feature data can be operative to provide feature data that expresses information for a plurality of segmented patient scans that have been normalized, aggregated, and compressed. The comparison module can be further operative to provide a confidence level for the categorization indicator. The system can further include a follow-up module operative to identify trends in changes to the categorization indicator over time. The system can further include an outgoing communication interface responsive to the comparison module and operative to provide the categorization indicator to a remote location. The source of patient imaging data can comprise part of an incoming communication interface. The system can further include an aggregate result analysis module operative to perform statistical analysis of results from the comparison module for a plurality of patients. The aggregate result analysis module can include correlative logic operative to determine relationships between treatment methods and categorization indicators for the plurality of patients. The system can further include a classification adjustment module responsive to the comparison module and operative to adjust the categorization information based on results from the comparison module. The disease characteristic categorization information can include a plurality of subcategories for different patient groupings. The subcategories can be based on symptom groupings. The subcategories can be based on patient demographic groupings. The feature data can further include further patient characteristics extracted from additional sources. The further patient characteristics can include patient demographic information. The further patient characteristics can include patient symptom information. The source of feature data can be a source of imaging data resulting from statistical analysis of image data acquired from the joints affected by different diseases. The source of feature data can be a source of imaging data resulting from dispersion analysis of image data acquired from the joints affected by different diseases. The source of feature data can be a source of imaging data resulting from principal component analysis of image data acquired from the joints affected by different diseases. The source of feature data can include categorization information for a global significant cartilage loss category and for at least one localized significant cartilage loss category. The source of feature data can include categorization information for a global significant cartilage loss category, for a global cartilage loss tendency category, and for at least one localized significant cartilage loss category. The source of feature data can include categorization information for a global significant cartilage loss category, for a global cartilage loss tendency category, for a medial condyle significant loss and medial trochlea loss tendency category, and for a medial condyle significant loss, trochlea significant loss, and medial condyle loss category. The source of patient imaging data can be operative to provide a fully automatically segmented imaging data set. The source of patient imaging data can be operative to provide a imaging data set having a resolution that is significantly lower than a resolution for the acquisitions on which the feature data is based. The source of patient imaging data can include digital identifiers associated with imaging data for a particular patient. The digital identifiers can include patient identifiers, physician identifiers, and joint identifiers. The source of patient imaging data can include error correcting codes. The error correcting codes can include codes associated with imaging data for a particular patient. The source of patient imaging data can include a format identifier associated with imaging data for a particular patient. The system can further include a patient verification module. The source of patient imaging data can be a source of magnetic imaging data resulting from a magnetic resonance imaging acquisition from the joint of the patient.

[0009] In another general aspect, the invention features a musculoskeletal imaging method that includes obtaining a patient imaging data set resulting from an imaging acquisition from a joint of a patient, obtaining feature data resulting from imaging acquisitions from joints of different individuals affected by different diseases, wherein the feature data includes extracted disease characteristic categorization information for a plurality of disease categories, comparing the patient imaging data set with the feature data, and providing, based on results of the step of comparing, a categorization indicator for the patient imaging data set that indicates a correspondence between information in the patient imaging data set and the extracted disease categories for which there is extracted categorization information in the feature data.

[0010] In preferred embodiments, the method can further include the step of determining whether to indicate a treatment for the patient based on the categorization indicator. The step of determining can be operative to determine whether to indicate the use of a drug therapy. The step of determining can be operative to determine whether to indicate the use of a disease modifying osteoarthritis drug. The step of determining can be operative to determine whether to indicate the local administration of an anti-inflammatory agent. The step of determining can be operative to determine whether to indicate the use of an exercise therapy. The step of determining can be operative to determine whether to indicate surgery. The step of determining can be operative to determine whether to indicate a rehabilitation method. The method can further include the step of determining whether to admit the patient to a clinical trail based on the categorization indicator. The method can further include the step of deriving from the categorization indicator and from other categorization indicators for other patients information about a treatment performed to at least some of the patients. The method can further include the step of determining whether to indicate a sports training regimen based on the categorization indicator. The method can further include the step of comparing the categorization indicator with results from other diagnostic methods. The step of obtaining a patient imaging data set can include obtaining a magnetic resonance imaging data set resulting from a magnetic resonance imaging acquisition from the joint of the patient.

[0011] In a further general aspect, the invention features a musculoskeletal imaging system that includes means for obtaining a patient imaging data set resulting from an imaging acquisition from a joint of a patient, means for obtaining feature data resulting from imaging acquisitions from joints of different individuals affected by different diseases, wherein the feature data includes extracted disease characteristic categorization information for a plurality of disease categories, means for comparing the patient imaging data set with the feature data, and means for providing, based on results of the step of comparing, a categorization indicator for the patient imaging data set that indicates a correspondence between information in the patient imaging data set and the extracted disease categories for which there is extracted categorization information in the feature data.

[0012] Systems according to the invention are advantageous in that they can enhance the amount of information available from a patient image. Because patient image data is compared with feature data that is selectively extracted from imaging data for a number of individuals in a particular population, systems according to the invention can be made to differentiate between subtly different types of patterns of disease progression. As a result, it may be possible to obtain meaningful diagnostic information about a patient's joint from a single image, instead of waiting to measure changes over time.

[0013] Early availability of diagnostic information is important in the case of progressive joint diseases, such as OA, because it can allow for early treatment of these diseases, and thereby help to delay or avoid their later stages. Simple changes in exercise habits or weight loss, for example, may delay more severe measures, such as surgery, if they are taken at the earliest signs of disease. It may even be possible to detect changes in the joint before the patient exhibits any symptoms, and thereby delay the onset of symptoms altogether.

[0014] Diagnostic information available from systems according to the invention may also be used to more precisely target treatments prescribed for particular individuals. Some types of patterns of progression may indicate that a patient is a good candidate for a particular type of drug treatment, for example, while others may show that such treatments do not warrant the risk of side-effects or that another treatment is most likely to produce the best results. This targeting information should result in earlier, more appropriate treatment for each of the different types of progression, and should help to reduce costs associated with treatments that are inappropriately targeted.

[0015] Systems according to the invention can also be highly objective. Because they can be made to rely exclusively on an automated comparison of image information with extracted feature information for a population, bias and human error can be kept to a minimum. This objectivity is very important in designing and monitoring scientific studies and clinical trials. As a result, it may be possible to more quickly and effectively develop additional treatments for joint diseases.

[0016] Systems according to the invention may also allow diagnosis to take place from lower resolution images. Because a comparison takes place with extracted feature information from a number of different patients, it may detect information too subtle to be seen by a radiologist. Less expensive, lower resolution MRI systems may therefore be used in some diagnostic and screening settings. Alternatively, shorter acquisition times can be used for a given imaging instrument, making the acquisition more comfortable for the patient and improving the availability and throughput of the instrument.

[0033] Referring to FIG. 1, an overall joint evaluation system 10 according to the invention can include one or more of the following modules: a patient interface 12, a recordkeeping module 14, a feature derivation module 16, a standard-setting module 18, a patient data processing module 20, and a diagnosis/treatment module 22. As will be discussed in more detail below, the modules and parts thereof can be implemented by different individuals and/or organizations. One of ordinary skill in the art would also recognize that it is also possible to combine, redistribute, or even omit functions for the various modules to achieve a different breakdown.

[0034] For example, recordkeeping functions can be performed by a larger organization, such as a teaching hospital, which also performs all of the other functions of the system. In smaller office settings, the recordkeeping might be performed by an outside service provider, or an affiliated hospital. And in some implementations, such as in a mobile screening clinic, it might even be possible to perform screening functions without performing any significant storage of patient data.

[0035] The patient interface 12 can include a patient exam data acquisition interface 24 and a patient image data acquisition interface 26. The patient exam data interface is designed to receive patient data, such as would be obtained during an exam by a physician. These data can include patient demographic data, such as age, race, occupation, height, and/or weight. The patient exam data can also include symptom data, such as self-reported patient data, results of physician-administered diagnostic tests, and/or other data reflective of the patient's symptoms. These data can be entered into the system directly, or it can be relayed from other systems.

[0036] The patient image data acquisition interface 26 receives image data reflective of spatial features of the joint of the patient, such as MRI data. These data can be obtained in a variety of ways, such as by being relayed directly from an imaging system, by being transferred electronically from another hospital, from a disk, or even by digitization of hard copies. Preferably, the evaluation is performed on a patient's knee joint, but other joints could also be evaluated.

[0037] The recordkeeping interface can include patient data storage 28 and feature data storage 30. The patient data storage can store some or all of the patient exam data and/or patient image data. The feature data storage can store feature data extracted from magnetic resonance imaging data resulting from magnetic resonance imaging acquisitions from joints of different individuals in a particular population of individuals. As will be discussed in more detail below, the feature data can include disease characteristic categorization information that allows a patient's image data to be categorized.

[0038] The patient processing module 20 can include a patient data processing module 40, a comparison module 41, a follow-up module 42, and a presentation module 44. The patient data processing module 41 performs operations to convert patient image data into a format that can be compared with the feature data. Its functions can include segmentation, sub-pixel interpolation, and data projection, although other types of operations may also yield suitable results.

[0039] The comparison module 41 performs comparisons between the processed patient image data and the feature data in order to evaluate the state of the patient's joint. In the present embodiment, the functions of the patient processing module and the comparison module are kept separate, allowing results of the patient data processing module to be stored. Although this functional separation is not necessary, storing processed image data allows for further comparison operations to take place without repeating patient processing operations. It may also allow the actual patient image data, which tends to be voluminous, to be stored offline.

[0040] The follow-up module 42 can combine information from the comparison module/or the patient data processing module with prior comparison information to obtain information about how the patient's joint has changed over time. The presentation module receives information from the comparison module 41 and/or the follow-up module 42 and formats it for review by a physician, the patient, and/or other subsystems, such as the data gathering module 32 and/or the diagnosis/treatment module 48. The presentation module can create result sets that take a variety of different forms, such as an enhanced image, one or more numerical values, or a suggested diagnosis and certainty value.

[0041] The patient diagnosis/treatment module 22 can include a diagnosis module 46 and a treatment module 48. The diagnosis module can provide a diagnosis indication and a corresponding confidence level for that diagnosis indication based on the information derived by the patient data processing module 20. The treatment module can provide a treatment indication and a corresponding confidence level for that treatment indication. Both the diagnosis indication and the treatment indication should preferably be in the form of a recommendation or a cross-check for the physician. It may be possible, however, to use these indications directly in certain circumstances. Truly blind clinical trial decisions could be made in an automated fashion based on information from the diagnosis module, for example, or automated rehabilitation equipment could derive operating parameters directly from the treatment module.

[0042] The feature derivation module 16 can include a data gathering module 32 and a data processing module 34. The data gathering module can include a set of functions designed to select and/or classify population data, such as MRI data. The data processing module 34 includes logic operative to derive feature data from the population data assembled by the data gathering module.

[0043] The standard setting module 18 can include a feature data standard setting module 36 and a data interpretation standard setting module 38. The feature data standard setting module sets standards for the acquisition of data and the processing of those data into standard sets. The data interpretation standard setting module sets standards for automated diagnosis and treatment operations.

[0044] As presented above, functions of the various system modules can be performed by different organizations, which may or may not be physically close. In one embodiment, one organization publishes one or more feature data sets, making them accessible to a number of evaluation organizations, such as hospitals or other health care providers. A third organization issues standards to ensure that the feature data sets are meaningful and actions taken based on them are sufficiently safe. Other organizations can provide further services, such as imaging services, recordkeeping services, or patient data processing services.

[0045] Some of these services can be provided over communications channels, such as the internet, through well-known interfaces such as Hypertext Transfer Protocol (HTTP) or extensible Markup Language (XML). A patient data processor, for example, could receive images from a client hospital through an secure internet connection, evaluate the images based on one or more feature data sets, and then transfer a result set back to the hospital through the same channel.

[0046] Referring to FIG. 2, overall operation of the joint evaluation system 10 begins with the acquisition of joint image data for a population (step 50). This joint data includes at least one image for each member of the population, although it may also include one or more follow-up images for these individuals, as well as demographic and symptom data. Members of the population are selected based on a desired set of characteristics. It may be desirable, for example, to sample the overall population of the region or country in which the system is being used, including both healthy patients and patients exhibiting ailments in their joints. It may also be desirable to replicate the characteristics of certain age groups. And particular specialists may also find it useful to work with a population made up of patients exhibiting certain types of symptoms or diagnoses.

[0047] The system then derives one or more feature data sets from the data acquired for the population (step 52). This derivation process extracts a subset of desired information from the patient data sets. The desired information would typically be that information which is most relevant to disease categorization, although other objectives might also be satisfied with different types of extraction processes. A number of different methods could be employed to perform the feature extraction, such as principal component analysis or dispersion analysis.

[0048] Once the feature data set has been derived and evaluated, it can be made available to one or more evaluation systems (step 54). The feature data sets may also be adjusted periodically or continuously (step 56). This adjustment can include continuing to acquire and process data from new or existing members of the population of interest, and/or it can take place through the incorporation of information from evaluation of patients based on earlier models (step 58).

[0049] Joint evaluation for a particular patient begins with the acquisition of joint imaging data (step 60). This acquisition can include the acquisition of images of one or both of the patient's knee joints using an MRI imaging system. The organization evaluating the patient can perform the acquisition, or it can receive the images from another source. MRI imaging parameters suitable for joint evaluation are described, for example, in U.S. Pat. No. 6,560,476, entitled “Evaluating Disease Progression Using Magnetic Resonance Imaging,” issued May 6, 2003, and herein incorporated by reference.

[0050] This methodology can be supplemented with a technique known as Active Shape Modeling (ASM). ASM is used to build a statistical shape model of the surface (or volume) to be segmented using semi-automatically segmented data as a training set. The model is composed of a mean shape and a shape variation space spawning the training set. Automatic segmentation is carried out, in the 3D image, by an alternating sequence of rigid and gradient based deformations applied to the mean shape. Before being applied to the segmented shape, deformations are projected into the model variation space so that resultant shape remains consistent with the training set. ASM is described in more detail in, for example. Active shape models—their training and application, by T. F. Cootes, C. J. Taylor, D. H. Cooper, and J. Graham, Computer Vision and Image Understanding, 61(1):389, January 1995, which is herein incorporated by reference.

[0051] The system then processes the patient data, including the joint imaging data (step 62). The system then compares the resulting processed patient data with the feature data (step 63). The result of this comparison operation can be a categorization indicator and a confidence value. The categorization indicator is an indication of which category the patient's joint best matches, and the confidence value is an indication of the level of confidence exhibited by the match. Also available is the underlying image comparison data.

[0052] If earlier joint data is also available, the system can derive trend information from the series of data sets (step 64). This derivation step can be part of an enhanced comparison operation in which the different joint data are processed in a single operation, or it can be a separate operation. The resulting data are then presented in a patient report (step 66), which can include a physician-readable portion, a patient-readable portion, and/or a machine-readable portion.

[0053] Based on the report, the physician, patient, and/or the system can take or refrain from taking action (step 68). Actions can include a wide variety of treatment-oriented actions including prescribing or adjusting a prescription of medication, surgery, exercise, or prosthetics. These types of actions are typically initiated by a physician or other health care practitioner, but in some cases, they may be automated. For example, a rehabilitation-oriented exercise machine may be able to adjust its resistance settings based on a patient report accessed electronically. And an implanted drug delivery system could also adjust drug delivery levels based on an electronically accessed patient report.

[0054] Other types of actions can also be taken, such as choosing to admit a patient in a clinical trial, certifying results for a clinical trial, transferring results to a database for further study, or comparing results to other types of diagnostic methods. In veterinary applications, actions might include choosing to breed an animal, selecting an animal for a particular task or training regiment, or even euthanizing an animal that cannot be treated. These types of actions will typically be performed by a health care practitioner, although they may also be automated or performed by less skilled personnel.

[0055] Referring to FIG. 3, the system 10 can apply different analytical techniques to different patient populations P1, P2, . . . PX. Each of these populations may exhibit the same or different categories C1, C2, . . . CN, . . . CO. The populations may be selected based on a variety of different criteria, such as age, occupation/activity level, or a preexisting diagnosis. The categories may be derived based on particular objectives, such as the selective diagnosis of different types of joint diseases. One population of interest is that of highly trained athletes, because their intensive training regimens may result in patterns of wear that differ from the general population. The health of their joints may also benefit from more frequent screening that do members of the general population.

[0056] Referring to FIG. 4, the patient data storage preferably includes patient record storage 80, and patient validation logic 82. The patient record storage can include storage for patient identification information 90, such as name, address, and identification number (e.g., social security number). The patient record storage can also include storage for patient history information 92, such as prior symptom and treatment information, prior patient image data, and prior processed patient image data. The patient record storage can further include storage for patient demographic information, such as patient age, race, and occupation. The patient record storage can also include storage for patient imaging data 96, such as one or more active MRI image data files. Error correction/detection encoding 98 may be provided at the record level, or it can be provided in the underlying systems.

[0057] The patient validation logic 82 helps to avoid errors in processing patient data. To this end it checks different parts of the patient records before allowing processing on a record to ensure that the correct record or the correct part of the record is being processed. To avoid comparing the knees of different patients with the same name, for example, the validation logic can flag situations where patient names match, but identification numbers, ages, or other information does not match. Similarly, to avoid comparing one patient's left knee with his or her right knee, a knee identification field can be compared with a simple test of the image data. Although the patient validation logic 82 is particularly important in the case of some ongoing clinical applications of the system, it may not be important or even required in other applications, such as one-time screening or veterinary applications.

[0058] The following example illustrates an approach to designing a system intended to allow for the assessment and diagnosis of osteoarthrisis (OA) using MRI cartilage layer thickness measurements. Although diagnosis based on cartilage thickness progression could still be made more precise, promising results have been obtained using the concept of normal joint templates. This approach employs a normalizing procedure that allows for the comparison of patient thickness maps to normal templates. Scaled difference maps can then be generated depicting cartilage deficit with respect to the normal reference.

[0059] To enable diagnosis based on normalized thickness for a particular patient, a decision rule is defined on the basis of synthetic quantitative parameters. One possibility is to compute the fraction of articular surface for which the scaled deficit exceeds a given threshold. The pathology severity could be related to this value and classes could be defined using different fraction thresholds. The main drawback of this approach is likely to be in the loss of topographic information, which may be critical, in the diagnosis decision rule.

[0060] The system design objective in this example is to create an objective, automatic and robust classification method—based on scaled difference—that takes into account spatial information. The general approach is to identify pathology patterns classes in an OA population, and this approach may be of great interest as a patient selection step for clinical studies. The identification of a few representative cartilage wear patterns may also constitute a solid base for pathology understanding, particularly if correlated with other indicators such as joint kinetic moments or meniscus tear scoring.

[0061] Signal classification is a known technique in “hard science” fields such as artificial intelligence, pattern recognition, computer vision, or statistical analysis. It can be decomposed in three main tasks: feature extraction, feature-based class definition and finally feature classification. In the present case, given that classes are not defined a priori, the problem is known as an unsupervised classification problem.

[0062] Feature definition has a strong influence on class definition and final classification. Because no trusted and explicit (parametrical, geometrical or anatomical) criterion were available for feature selection, it was decided to use a general feature definition based on information theory and statistical analysis. A basic framework consists of following steps:

image information compression, keeping meaningful data structures using methods such as principal component analysis (PCA), factor analysis (FA), or independent component analysis (ICA)

• • [0063] image classification in the reduced space (with supervised or unsupervised methods, that may or may not use underlying probabilistic models)

[0064] A simple classification procedure based on PCA and K means algorithm for clustering was tested. It is noted that this approach does not rely on an explicit underlying data model, but such a model could be derived through statistical non-linear mixture modeling.

[0065] Femoral cartilage thickness maps of 110 OA patients and 10 normal age-matched patients were first processed from cartilage interface contours segmented with a semi-automatic method, as presented in U.S. Pat. No. 6,560,476, discussed above. The maps employed a cylindrical parameter space using a fitted cylinder as a primitive. A normal thickness template was first generated by aligning, scaling and averaging the maps for the 10 non-arthritic subjects (see FIGS. 9-10).

[0066] Scaled difference maps could then be generated by subtracting patient OA maps from the normal template and scaling the resulting difference by a standard deviation (SD) template. Referring to FIG. 6, for example, a map of mean scaled differences among members of an OA population of 110 subjects shows average deficit and surplus cartilage thickness in a the population. Medial condyle (left) and trochlea regions show clearly a significant deficit (>=3 SD) whereas other regions are not significantly thicker than normal.

[0067] To simplify computations, patient surface scaling and alignment were applied directly in the parameter space. Five different alignment and scaling methods have been evaluated, and their main characteristics are summarized in Table 1. The first four methods share the same principle: alignment of map bounding box centers, surface matching by uniform and linear scaling on the two parameter dimensions, bicubic thickness interpolation, and, finally, thickness scaling using the same factor as the one used for surface matching. The methods differ in the way the scaling factor is estimated: map anterio-posterio hight (AP), map medial-lateral width (ML), fitted cylinder radius (R), and cartilage surface square-root (SSR).

[0068]

TABLE 1
Summary of alignment and scaling methods
			Thickness	Thickness
	Alignment	Surface scaling	interpolation	scaling
Method 1	Alignment of	Uniform scaling along	Bicubic	Normalization
	the maps	the two parameter	interpolation in	of thickness
	bounding box	dimensions using the	the parameter	values by the
	centers	anterior-posterior	space	scaling factor.
		height of the map as
		the scaling factor (AP)
Method 2		Same procedure using
		the map medial-lateral
		width as the scale
		factor (ML)
Method 3		Same procedure using
		the fitted cylinder
		radius as the scale
		factor (R).
Method 4		Same procedure using
		the surface square root
		as the scale factor
		(SSR).
Method 5	Biquadratic surface warping based on		Surface
	feature points extracted from map		square root
	contours.		scaling factor
			(SAF)

[0069] The last method is a warping method based on feature points extracted from the map contours. Four anatomical points P1-P4 are first defined on the femur map as shown in FIG. 5. P1 and P2 are the local maxima points on the condyle parts, and P3 is the local minimum in between. P4 is the intersection between the bottom contour and the Blumensat line.

[0070] Other feature points are defined by equal subdvision of the contour sections defined by the anatomical features. Averaging these features among all the subjects defines the reference features. A biquadratic transformation is estimated for each subject to transform subject features to reference features. This transformation is then used to process the matched thickness maps. Finally, as for former methods, a scaling factor (equal to the surface square root ratio) is applied uniformly to thickness values.

[0071] The classification method provides an objective and topographical base for the pathology map class definition. Given the high dimension of data (each map contains typically 50000 points) and the great variety of scaled difference map topographies, information should be simplified and summarized (compressed or coded) to extract only the most significant patterns. Hence, feature definition and extraction is the first step for map (i.e., signal) classification. There is no explicit and objective formulation of these features because pathology classes are not defined. Yet, it seems likely that classes will rely on wear region pattern characteristics, including the number of regions, their positions, and their topographic properties.

[0072] One way to bypass explicit feature definition is to consider maps as signals and features as the uncorrelated information they contain. From this point of view, feature extraction is seen as a signal quantization process.

[0073] In this example, it was decided to follow a classical framework used in multivariate data analysis and information theory. Each map is viewed as a multi-dimensional set of variables. Information is first reduced exploiting linear correlation via principal component analysis (PCA), and then classes are defined using an unsupervised classification method applied on uncorrelated coefficients. PCA is actually an appropriate tool to code and decorrelate difference map information. It will naturally associate correlated zones and constitute an optimal base to automatically define variation classes. Class definition using PCA is optimal in the sense that classes best represent (relative to a chosen inter-map distance measure) the variance of the reference maps.

[0074] PCA application to image analysis has been successfully used in the computer vision field for face recognition. In this application, maps are considered as 2D images in a M×N dimension space. Each image of this space can be represented by column vector of length d=M.N, which is simply the lexicographic reading of that image. The result is an d-dimensional space, called map space, in which wear patterns can be characterized.

[0075] Referring to FIG. 7, a typical set of thickness scaled difference maps can be chosen as a reference to describe the variety of pathology cases. The PCA procedure defines a sub-space that contains the most significant information (i.e., variations around an average map) from the reference maps, automatically associating correlated map points, in order to build few typical base maps (eigenmaps) mostly representative of the main variations of the reference maps. Each reference map can then be represented, in a simplified way, as a linear combination (eigenmap components) of few eigenmaps. Only the main information is retained.

[0076] In mathematical terms, an eigenspace is spaned by an orthogonal base of eigenmaps defined as the eigenvectors of the correlation matrix of the reference signals. Eigenmap components are defined as the projection of the reference maps onto the eigenmaps. Eigenvectors are built orthogonally one to each other, which means eigenmap components are uncorrelated. They also maximise the variance of eigenmap components. Each eigenmap is associated to an eigenvalue proportionnal to the information quantity it contains.

[0077] Let T be the (d×q) matrix containing q reference maps represented as (d×1) column vectors t_i. In fact, PCA provides a geometrical solution to approximate each t_i by its orthogonal projection ({circumflex over (t)}_i) in a sub-space. This sub-space, spanned by (d×q) orthogonal matrix A, is chosen to maximize variance of ti components in it. In order to satisfy maximal variance property, A is composed of the q most significant eigenvectors (eigenmaps) of the correlation matrix C of T (C=(T− T).(T− T)^T). Let μ be the mean reference maps. {circumflex over (t)}_i can then be written as,
{circumflex over (t)}_i=AA^T(t_i−μ)+μ
Defining (q×1) vectors s_i (eigenmap components) as the components of the maps in the eigenspace, we have,
s_i=A^T(t_i−μ) and {circumflex over (t)}_i=As_i+μ
Moreover, defining the reconstruction error vector ε_i as:
ε_i=t_i={circumflex over (t)}_i
In other words, PCA provides a linear model for original data:
t_i=As_i+μ+ε_i
Data can even be more compressed by selecting only the q′ most significant eigenvectors. (q×q′) matrix A′_i, with q′<q, replaces matrix A.

[0078] Given the nature of wear process, data points have been observed to be highly spatially correlated. Most of maps show large regions where the signal can be approximated by a linear model. In other words, the signal is not highly textured and does not appear to need higher-order moment data models. Moreover, as far as diagnostic and classification are concerned, constant mean region patterns are sought and it is expected that PCA will simplify information and highlight these patterns. Many non-linear PCA extensions have been designed (probabilistic PCA, ICA . . . ) for other applications but it is not clear that they are necessary the present type of signal and issue.

[0079] Eigenmaps represent population variations along an associated axis of the average map. They are shown in FIG. 7 in decreasing order of variation importance. To each eigenmap is associated an eigenvalue that is equal to the standard deviation of the map population projected on the variation axis. Eigenvalues are propotionnal to the quantity of information the eigenmap can represent. Their “representativeness” is expressed in percentage points in FIG. 7. FIG. 8 shows information distribution among the thirty main components. The right plot represents the cumulative information they retain. Hence it is noted that the ten first components represent more then 80% of the information.

[0080] The first eigenmaps can describe great variations (large zones and/or high variation amplitude). The latest ones allow the system to model smaller details (small zones and/or low amplitude variations). The observable patterns or structures (groups of connex points having roughly the same value) represent thinner (or thicker) zones that allow to describe the most different maps. Given the nature of the map population signal, eigenmaps show simple and regular patterns like homegenous regions. Note that region sizes decrease with the component representativeness. In other words, selecting a restricted number of components allows to generate maps with a certain “spatial resolution”. The more components used, the smaller is the spatial resolution. This property will be considered for selecting the number of classification components.

[0081] Contrary to supervised classification cases, class definition is not treated as available in this case. Some clinical features, such as meniscus tear score or pain test results, could be used for classification, but the goal of the method is actually to propose a way to define classes on the base of quantitative cartilage thickness information.

[0082] The general concept of class definition raises the following questions:

• • [0083] How to decide when data variations are intra-class variability or inter-class variability
  • [0084] Variability definition itself raises the question of how to define distances between individuals and more generally how to model data
    In this case, data variations are caused by both spatial structure variability (point correlation variations) and point intensity variability. A map class modeling was defined in order to adequately describe these variations in terms of intra-class and inter-class variability.

[0085] The PCA compression step already assesses an implicit constraint on class definition because orthogonal projection eradicates the smallest data variations, in the sense of Euclidian distance between signals. The choice of the projection does not guarantee that compression is optimal in terms of classification performance. In contrast, methods using an a priori class definition, on a training set, can explicitly optimize projection axes so that inter-class variability is maximized (discriminant analysis).

[0086] The K-means algorithm was initially chosen in the map space because of its simplicity and efficiency. K-means is simple and powerful clustering algorithm that minimizes the average intra-class distance also called quantization error E. In other words, for a set of N eigenmap component vectors s_i, K-mean defines M class centers ω_k that minimize E(ω):

[0087] $E\left(\omega \right)=\sum _i\frac{1}{2}\underset{k}{\min }\left[\mathrm{dist}\left(s_i,{\omega }_k\right)\right],$
where dist(s_i, ω^k)=(s_i−ω_k)² in the Euclidian case.
The choice of Euclidian distance implies that classes have hyper-spheroid shapes. Other distance types like Malhanobis or generalized distance also exist allowing for seeking of hyper-ellipsoid and probabilistic shapes. Given the absence of specific knowledge about the class shapes, the Euclidian distance was chosen because of stability results. A combination of batch and on-line versions of the algorithm were used in order to optimize convergence speed and accuracy.

[0088] Once classes have been defined on the training set, a new map can be simply classified by projecting it into the map space and labeling it with class that is the closest to it, using Euclidian distance.

[0089] The method needs mainly the selection of two parameters: the number of significant components and the number of classes. Many different methods can be found in the litterature to select the optimal values, depending on the data model (Gaussian mixture, geometrical . . . ) and optimization criterion (MAP, ML, ratio intra-inter class distance . . . ). Some methods using comparison with simulated reference populations (GAP statistics) have also recently been proposed, but experimental results did not appear to be well adapted to the problem at hand.

[0090] Practically, the first parameter is not critical and can be heuristicaly assimilated to the detail level at which maps are coded before clustering step. The more components are selected, the smaller and less correlated will be the reconstruction error. On the other hand if the detail level is too high, the clustering step stability will decrease and class definition will be less significant. Initially the choice to use the 10 first principal components was made heuristically on the base of information quantity (@ 80% of total variance) and pattern detail level presumed to be relevant for pathology class definition. Beginning with a small number of components then allows for a preliminary low resolution classification that describes large-scale tendencies in a robust way.

[0091] The second parameter is more complex to optimize. The main pragmatic criteria were stability and robustness towards segmentation (human factor) and clustering (influence of initialization step) variabilities. A relatively homogenous mixture proportion is also important. Clusters with very few members are undesirable. Given these considerations five clusters were used for final classification results of the example population.

[0092] A normal template has been generated using 10 non-arthritic subjects. This population is relatively small because it was found to be important to select subjects in order to match normal population and OA population ages, but it should be straightforward to perform repeat studies for normal templates based on larger populations.

[0093] Refering to FIGS. 9 and 10, template thickness mean (FIG. 9) is relatively independent of the scaling method, but this is not the case for the standard deviation (SD) of mean surface and thickness map (FIG. 10). Sum of root-mean-square of standard deviation (SD-rms) is used as an indicator of scaling factor adequacy (see Table 2 and Table 3). It was noted in preliminary studies that SD-rms is very much dependant on the number of maps used to compute a template. Because anterior-posterior scaling factor (AP) gave relatively good results for different populations and because it was refered as a pertinent scaling method in the literature, it was selected for use in the following steps.

[0094]

	TABLE 2
	Method 1	Method 2	Method 3	Method 4
	(AP)	(ML)	(R)	(SSR)
Thickness SD	2.7240	2.8689	2.7668	2.5516
rms
(mm)
Thickness CV	11.0252	11.4131	10.9338	10.3995
rms
(%)

[0095] Table 2 presents root-mean-square of cartilage offset standard deviation and coefficient of variation for 10 normal age-matched femurs using four different scaling methods. In the table, AP stands for anterior-posterior, ML stands for medial-lateral, R stands for radius of fitted cylinder, and SSR stands for surface square root.

[0096]

	TABLE 3
	Method 1	Method 2	Method 3	Method 4
	(AP)	(ML)	(R)	(SSR)
Thickness SD	0.4151	0.4009	0.4400	0.4071
rms
(mm)
Thickness CV	24.4235	23.7568	25.5213	23.8798
rms
(%)

[0097] Table 3 presents root-mean-square of cartilage thickness standard deviation and coefficient of variation for 10 normal age-matched femurs using four different scaling methods. In the table, AP stands for anterior-posterior, ML stands for medial-lateral, R stands for radius of fitted cylinder, and SSR stands for surface square root.

[0098] Eigenmaps have been computed using baseline acquisition of 110 OA subjects normalized by a template based on 10 non-arthritic subjects. All maps have been projected in the eigenspace and reconstructed in the full dimension space. As an indication, FIGS. 13 and 14 show results for subjects 1 to 5 (columns, FIG. 13) and subjects 6 to 10 (columns, FIG. 14). Rows show reference map (scaled difference), projection in the eigenspace using 10 components, reconstruction error and associated class map (6 classes).

[0099] Projected maps are very similar to the original ones. The error information for each subject is relatively uncorrelated, which means that 10 components are sufficient to capture original signal structures relatively well. Finally it is noted visually that original global map structures are very similar to the ones present in the class map they have been assigned to. Classes therefore visually represent the main features of the signal well.

[0100] Three type of results have been generated:

• • [0101] baseline and follow-up classification comparaison
  • [0102] meniscus tear scoring and classification dependancy analysis
  • [0103] subjet gender and and classification dependancy analysis

[0104] Baseline and follow-up classifications have been carried out as a preliminary test to evaluate the consistency over time of classification changes. It is expected that the method will filter segmentation noise and provide a more significant general cartilage change indicator. Dependency between clinical features like meniscus tear score or subject gender allows the evaluation of class labeling consistency from a pathological point of view.

[0105] FIG. 15 presents scaled baseline difference maps of the five classes clustered by the K-means algorithm applied on 105 OA patients for which follow-up was available. Class a posteriori probabilities of the classes were computed. These class maps could qualitatively be interpreted as follows:

• • [0106] Class 1: global significant cartilage loss
  • [0107] Class 2: global cartilage loss tendancy
  • [0108] Class 3: medial condyle significant loss, medial trochlea loss tendancy
  • [0109] Class 4: medial condyle significant loss, trochlea significant loss and medial condyle loss tendancy
  • [0110] Class 5: very similar to class 3

[0111] Table 4 and Table 5 present scaled difference map classification probability evolution from base line to follow-up. Follow-up scaled different maps are classified in the eigenspace by minimizing Euclidian distances to the baseline class centroids.

[0112] Table 5 shows the class probability for baseline and follow-up. The main change seen is that patients for class C2 moved to class C1 which corresponds to a general cartilage loss. Table 4 presents classification transitions and their occurrences. Corresponding cartilage evolution is also listed. Note that 96 out of 105 are stationary, eight transitions correspond to a cartilage loss and only one reflects cartilage gain.

[0113]

	TABLE 4
	P(C1)	P(C2)	P(C3)	P(C4)	P(C5)
Base Line	0.3619	0.3143	0.1619	0.1238	0.0381
Follow-up	0.4095	0.2762	0.1524	0.1333	0.0286

[0114]

TABLE 5
Base Line	Follow-up
Classification	Classification	Cartilage evolution	Occurences
Classe i	Classe I	stationary	96
Classe 2	Classe 1	general loss	5
Classe 3	Classe 2	light lateral loss	1
Classe 3	Classe 4	medial and trochlea	1
		loss
Classe 4	Classe 3	medial and trochlea	1
		gain
Classe 5	Classe 4	medial and trochlea	1
		loss
			= 105

[0115] Meniscus tear scoring was carried out for each patient. Using a combination of anterior/posterior (anterior, medium, and posterior positions) and lateral/medial (lateral and medial position) axes, six meniscus regions were defined (lateral region codes: AL, ML, PL; medial region codes: AM, MM, PM). Each region was scored by a binary tear score (1=tear, 0=no tear). In order to simplify and filter information, data was summed up in four major regions (medial, lateral, anterior, posterior coded M, L, A and P—see FIGS. 16A and 16B). Major region scores were computed heuristically as follows in order to capture possible lateral or longitudinal gait tilt tendencies:

• [0116] M=1 if AM+MM+PM>1, else M=0
• [0117] L=1 if AL+ML+PL>1, else L=0
• [0118] A=1 if AM+AL>1, else A=0
• [0119] P=1 if PM+PL>1, else P=0

[0120] Tear conditional probabilities of major regions were evaluated a posteriori for the four more significant classes (see Table 6). Probabilities were analyzed in order to detect rough tilt tendencies. Results are therefore compared by pairs to evaluate tilt along longitudinal (region A vs. region P) and lateral (region M, region L) axis. Intuitively, it would be roughly expected that a high tear conditional probability in a specific region relative to its paired region should have an impact on the class map topography. For example, high tear probability in the medial compartment (M region) associated with a lower probability in the lateral compartment for a given class should be associated with larger and deeper wear patterns in the medial region for a corresponding class map. Because wear processes follow a complex and multifactorial evolution, this semi-qualitative analysis has been used very carefully as a classification coherence indication.

[0121]

	TABLE 6
	Tilt axis:
	Lateral	Longitudinal
	Tear location:
	Medial (M)	Lateral (L)	Anterior (A)	Posterior (P)
P(tear)	0.5238	0.5429	0.3810	0.4000
P(tear/C1)	0.3421	0.6316	0.1842	0.4737
P(tear/C2)	0.5152	0.4545	0.4242	0.3636
P(tear/C3)	0.7647	0.5294	0.6471	0.2941
P(tear/C4)	0.6154	0.6154	0.5385	0.4615

Lateral (0.5238 vs. 0.5238) and longitudinal (0.3810 vs. 0.4) tear probabilities for the whole population are relatively balanced. Medial tear probability for class C1 (0.3421) is smaller than the lateral one in that region (0.6316). Along the other axis, posterior probability is higher than the anterior one. Class C2 is relatively balanced for both axes. Class C3 is much more unbalanced in both axes (higher medial and anterior tear probabilities). Finally class C4, like C2, is relatively balanced along both axis.

[0122] Gender conditional classification probabilities are presented in Table 7. It is know from epidemiology litterature that women are more susceptible to OA. Hence it should be expected that female distribution among classes would be more focused on damaged classes. A possible sorting of classes in increasing pathological gravity order could be: C2, C3, C4 and C1. Note that females are mostly concentrated in class C1 (highest gravity) whereas male are concentrated in class C2 (lowest gravity).

[0123]

	TABLE 7
	Classes:
	C1	C2	C3	C4	C5
P(Ci/Male)	0.2286	0.4286	0.2286	0.0571	0.0571
P(Ci/Female)	0.4265	0.2647	0.1176	0.1618	0.0294

[0124] The present invention has now been described in connection with a number of specific embodiments thereof. However, numerous modifications which are contemplated as falling within the scope of the present invention should now be apparent to those skilled in the art. For example, the techniques described may be used in veterinary applications or for the imaging of other types of structures in the body. Another approach to normalization of joint structures is to employ landmark extraction and dense mesh correspondence. Therefore, it is intended that the scope of the present invention be limited only by the scope of the claims appended hereto. In addition, the order of presentation of the claims should not be construed to limit the scope of any particular term in the claims.