An Atlas-based Approach to Study Morphological Differences in Human Femoral Cartilage between Subjects from Incidence and Progression Cohorts: Mri Data from Osteoarthritis Initiative

Recent advances in magnetic resonance (MR) imaging have made it the modality of choice for assessment of cartilage status for osteoarthritis (OA) (Roemer et al., 2009; Eckstein, 2007). MRI allows assessment of morphological changes using high resolution 3D imaging as well as evaluation of changes at the biochemical/molecular level using MR relaxometry techniques (T2: spin-spin relaxation and T1: spin-lattice relaxation with spin locking) and delayed Gadolinium enhanced MRI of cartilage (dGEMRIC) (Blumenkrantz & Majumdar, 2007). Research on the potential application of MRI to detect and classify osteoarthritis and monitor progression is an area of intense research (Hunter et al., 2009). A number of review articles detail the technical limits, accuracy and precision of MRI including the choice of pulse sequences and image analysis methods for morphological assessment in osteoarthritis (Raynauld, 2003; Eckstein et al., 2009; Xing, 2011). Cartilage morphology metrics from MRI are now being actively explored as imaging biomarkers of disease status and progression. The sequences that provide the best contrast of cartilage to adjacent tissue as well as the best resolution are the fat-suppressed spoiled gradient echo and fast double echo and steady state (DESS) with water excitation (Blumenkrantz & Majumdar, 2007; Roemer et al., 2010). Segmentation of the cartilage is a challenging task and few completely automated segmentation algorithms have been proposed. A fully automated segmentation algorithm based on Active Shape Modeling (Fripp et al., 2007) has been implemented and evaluated for normal cartilage. A more recent automated algorithm that is extensible to osteoarthritic cartilage first segments the bone-cartilage interface and then classifies cartilage voxels based on texture analysis. However, the large-scale applicability of these algorithms to normal and diseased cartilage is not yet demonstrated. Semiautomated algorithms using region growing, edge detection, and spline fitting have been used successfully in large-scale projects (Raynauld, 2003; Duryea et al., 2007). Morphological indices used to classify disease status and progression, include global and regional measures of total cartilage volume (raw and normalized), cartilage surface area, cartilage thickness (global and regional), and denuded area (Pelletier et al., 2007; Eckstein et

detect biochemical changes will be essential to the early diagnosis and treatment of pathologies.MRI affords this ability through measurement of MR indices that are sensitive to local biochemical changes.The MR indices that have been established as sensitive markers of osteoarthritis are: T2, T1 and dGEMRIC (Blumenkrantz et al., 2007).Findings from these studies have shown these parameters to correlate with clinical findings and a potential to uses these parameters in a clinical environment to detect and monitor OA.These techniques have the power of detecting cartilage loss at an early stage when the damage is still reversible (David-Vaudey et al., 2004;Li et al., 2005;Phan et al., 2006;Regatte et al., 2006;Koff et al., 2007;Li et al., 2007;Link et al., 2007).T2 relaxation time (spin-spin relaxation) is a noninvasive marker of cartilage integrity as it is sensitive to tissue hydration and the collagen-proteoglycan matrix (Blumenkrantz et al., 2007;Taylor et al., 2009;Link et al., 2010).In normal cartilage, T2 is low due to the immobilization of the water protons in the collagen-proteoglycan matrix.Depletion of collagen and proteoglycan in osteoarthritis renders protons more mobile which is reflected as longer T2 relaxation times.Another factor that contributes to the increased T2 is due to an increase in water content in diseased cartilage.However, the link between T2 and biochemical changes is complex; T2 changes have been correlated with changes in collagen content but not with proteoglycan loss (Regatte et al., 2002).An interesting feature of cartilage T2 is the spatial distribution within the cartilage with a steep gradient from the subchondral bone to the cartilage surface.Prior studies have explored the relationship between T2 and cartilage morphometry: T2 has been found to be inversely correlated with cartilage volume and thickness (Dunn et al., 2004).Studies analyzing MRI T2 relaxation time have also been performed on the data from OAI (Carballido-Gamio et al., 2011;Pan et al., 2011).Carballido-Gamio et al. propose a new improved technique that involves flattening of the cartilage and application of texture analysis to create T2 maps.The texture analysis performed using gray-level co-occurrence matrix in direction parallel and perpendicular to the cartilage layers.The longitudinal analysis is performed on data at baseline, 1-year follow-up, and 2-year follow-up.The results show several textures features showing higher values perpendicular to the cartilage layers where other texture feature exhibiting higher values parallel to the cartilage layers.These finding provide initial results to an alternative approach to study cartilage in OA (Carballido-Gamio et al., 2011).Another study investigates the correlation between T2 values obtained from MR images at 3.0T and the vastus lateralis/vastus medialis cross-sectional area (VL/VM CSA) in a group of middle aged subjects from the incidence cohort (non-symptomatic).The authors used axial T1W images to calculate the CSA values of the thigh muscles.The results showed that higher values of T2 correlated with greater loss of cartilage in OA.Regression values indicated that the VL/VM CSA values are inversely proportional to the T2 values (Pan et al., 2011).T1 refers to measurement of T1 in the rotating frame and probes very low frequency interactions such as between water and large macromolecules (like proteoglycan) that are present in the extracellular matrix of the cartilage (Blumenkrantz et al., 2007;Taylor et al., 2009).T1 has been shown to be correlated with proteoglycan disruption and shows regional variations similar to T2.Studies have shown that cartilage T1 values in OA subjects to be increased compared to controls that is in accordance with the model of proteoglycan disruption with disease.Comparing T2 and T1, Majumdar et al. found T1 to be more sensitive than T2 in distinguishing OA cartilage from healthy cartilage (Majumdar 2006).
In dGEMRIC studies, the charged contrast agent GdDTPA 2-is injected intravenously into a patient and allowed to diffuse into articular cartilage and imaging is performed 2-3 hours after contrast administration (Blumenkrantz et al., 2007;Taylor et al., 2009).The proteoglycan molecule consists of a central protein core to which a large number of negatively charged side chains called glycosaminoglycan (GAG) are attached.The GAG side chains repel the negatively charged contrast agent so that regions of low GAG concentration (as in diseased cartilage) show greater contrast uptake than normal cartilage.The distribution of the contrast agent (and GAG) is calculated from a T1 map of the cartilage.Like T2 maps, the T1 maps are heterogeneous indicating that the spatial distribution of GAG is heterogeneous.The dGEMRIC has been used to evaluate osteoarthritis in knee cartilage and T1 values correlate to the KL grading scale (Williams et al., 2005).dGEMRIC has shown potential in early detection of osteoarthritis as it provides direct information on the distribution and content of GAG in cartilage.From the discussion so far, it is clear that regional quantitative assessment of morphological/biochemical cartilage is more sensitive than global changes, and the analysis requires one to account for spatial heterogeneity in the indices as well as the variability of these indices in a normal population.To address these factors, studies have used T and Z scores successfully to study cartilage loss, T2, T1 and dGEMRIC changes in OA (Blumenkrantz et al., 2007).T and Z scores provide an estimation of the difference between patients and healthy young subjects and the difference between patients and their age matched healthy subjects respectively.Burgkart et al. showed in their study that the cartilage volume measurements when normalized on joint surface area before diseased state increase the accuracy and applicability of T and Z scores by reducing inter-subject variability (Burgkart et al., 2003).However, the T or Z score are still based on the mean value for normal subjects in a cartilage compartment.Thus, it averages the spatial distribution of the index over the compartment and given the spatial variation of the MR indices, this approach may mask subtle changes.Morphological quantification of osteoarthritis currently includes measurement of subregional cartilage loss.Detailed studies have revealed a small loss of 1-2% in cartilage thickness annually and a high degree of spatial heterogeneity for cartilage thickness changes in femorotibial sub-regions between subjects (Eckstein et al., 2010;Frobell et al., 2010;Wirth et al., 2010).Most quantitative studies divide the cartilage into a few compartments and report cartilage thickness (mean and standard deviation) for each compartment.Wirth et al. reported a technique for regional analysis of femorotibial cartilage thickness based on quantitative MRI (Wirth et al., 2008).The latter paper uses an elegant algorithm driven identification of sub-regions with user-controlled parameters to define the sub-regions.Wirth et al proposed that these parameters could be tuned according to the regional cartilage progression with OA.This technique represents a significant advancement in the automated and quantitative assessment of cartilage thickness.However, localized changes smaller than the size of the sub-regions may escape detection even with this technique.Eckstein et al explored the magnitude and regional distribution of differences in cartilage thickness and subchondral bone area associated with specific Osteoarthritis Research Society International (OARSI) JSN grades.Their regional analysis provided quantitative estimates of JSN related cartilage loss and revealed that the central part of the weightbearing femoral condyle as the most strongly affected (Eckstein et al., 2010).In a recent paper, Wirth et al extended their regional analysis to identify spatial patterns of cartilage loss in the medial femoral condyle.They localized the posterior aspect of the central, weight-bearing medial femoral condyle as showing the greatest sensitivity to change in disease status (Wirth et al., 2010).These regional analyses clearly indicate that cartilage loss occurs in a spatially heterogeneous manner and sometimes, in very localized regions.Stammberger et al. proposed elastic registration of 3D cartilage surfaces to detect local changes in cartilage thickness; in both synthetic and volunteer data, thickness differences recovered from the registration method were similar to that from using Euclidean distance transformations (Stammberger et al., 2000).Cohen et al. generated templates of cartilage of the patellofemoral joint and demonstrated the potential of using the standard thickness maps by comparing it with thickness maps generated for individual patients to identify regions with maximum loss of cartilage in patients with Osteoarthritis (Cohen et al., 2003).Recently, Carballido-Gamio et al. performed inter-subject comparison of knee cartilage thickness after registration to a common reference space (Carballido-Gamio et al., 2008).They measured the thickness at each point on the bone-cartilage interface and used affine and elastic registration techniques for point-wise comparison of cartilage thickness.They established the reproducibility of this method for intra-and inter-subject cartilage thickness comparisons and concluded that the techniques could be used to build mean femoral shapes and cartilage thickness maps.They showed that the proposed technique was an accurate and robust way to analyze inter-subject thickness point wise.However, there are no reports to date on the construction of a cartilage atlas to automatically characterize cartilage morphology in a population group (i.e., a cartilage atlas), which can be used to assess localized morphological or parametric differences between cohorts.Our approach extends prior work on creation of standard maps of cartilage thickness to an atlas-based approach.The atlas-based approach allows automatic detection of voxel based differences between cohorts (e.g., segregated by age, gender, ethnicity, disease status) as well as enables tracking of longitudinal changes.The approach is general and differences are not restricted to morphological differences at the voxel level but extend to parametric maps of T2, T1, and dGEMRIC.The approach has the potential to allow automatic and accurate detection of subtle changes as it leverages the statistical power of the cohort size.Thus, it is possible to model variability within the cohort and identify significant differences between the cohorts at the voxel level.This unique approach allows measuring localized changes in cartilage morphometry without any previous knowledge of probable regions of changes.The approach is also not restricted to detecting differences/changes in cartilage; it can also be extended to other anatomical structures that are impacted by osteoarthritis (e.g., bone).In this chapter, we outline the atlas based method and apply it (i) to detecting disease-based differences in cartilage morphometry, and (ii) to creating a bone atlas from knee MRI of normal healthy young adults.

Data selection criteria and MRI
Data for the current project was obtained from the Osteoarthritis Initiative.Realizing the importance of OA and its impact on millions of patients (socially and economically) National Institute of Health (NIH) started a four-year observational study called the Osteoarthritis Initiative (OAI).This initiative was lead by the National Institute of Health (NIH) in collaboration with an academic and industrial consortium to obtain clinical, biochemical and imaging data on a large cohort of subjects to follow the onset and progression of osteoarthritis.The OAI provides an extremely rich database of high resolutions MR images of human knees and is made publicly available to the research community (www.ucsf.edu).OAI provides MR images for 4,796 participants (ver.0.E.1) at baseline and the follow-up data for these participants at periodic intervals through other data release versions.Availability of this rich dataset, have resulted in several promising research methodologies and significant findings.This study uses sagittal knee magnetic resonance images made available in OAI data release version 0.A.1.A water-excitation double echo steady state (DESS) imaging sequence was used to generate the MR images on a Magnetom Trio, Siemens acquisition machine at 3.0T.The imaging parameters used in the acquisition protocol are repetition time (TR)/echo time (TE): 16.3/4.7ms,x and y resolution: 0.365, slice thickness: 0.7 mm and field of view (FOV): 140 mm.The OAI data version 0.A.1 for images (MRI and X-ray) refers to the baseline line data and represents its first public release.This release contains 200 incidence and progression cohort participants.The participants in this group were further stratified based on gender and clinic (four recruitment centers).OAI adopts the following criteria to classify subjects as belonging to the progression or incidence sub-cohort.The progression cohort consists of participants that have symptomatic OA at the onset of the study.The incidence cohort contains participants with no prior symptomatic OA at the start of the study, but is at a higher risk of knee OA.OAI defines symptomatic knee OA in participants if they meet both of following criteria: 1. Knee symptoms defined by "pain, aching, or stiffness in or around the knee on most days" for at least a month in the past year and 2. Radiographic evidence of tibiofemoral osteophytes that is comparable to Kellgren and Lawrence grade  2. Participants with no symptomatic knee OA (as defined above) in either knee at baseline, but exhibit characteristics such as frequent knee symptoms with no radiographic OA or meeting two or more eligibility factors, are classified in the incidence sub cohort.The cartilage atlas was created from 30 participants chosen from the available set of 200 incidence and progression cohorts.All atlas subjects belong to the incidence cohort and are Caucasian males.Further, subjects in the incidence group with KLG scores in the range of 0-2 ('absent' to 'mild' OA) were used to create the atlas.Sagittal knee MR images of the selected 30 male participants were used to create the femoral cartilage image atlas.The demographics of the participants are as follows: mean  standard deviation of the age 66.5  8.2 years and mean  standard deviation of the KLG grade 1.56  1.1 (more details can be obtained from our previous publication (Tameem et al., 2011).In addition to the 30 participants chosen to create the cartilage atlas, 10 male participants were chosen from the incidence and progression group each to compare localized changes in femoral cartilage with disease condition.It should be noted that the ten incidence cohort subjects chosen for the comparison were not part of the atlas cohort.The mean  standard deviation of age for participants chosen from the incidence and progression cohort was 64 ± 8.87 and 67.5 ± 8.78 respectively.The average knee osteophyte (KOST) grade, average knee joint space-medial (KJSM) grade, average knee joint space-lateral (KJSL), and Kellgren and Lawrence (K&L) grade for the incidence cohort is 1, 0, 0, and 1 respectively.The average KOST, KJSM, KJSL, and K&L grade for participants from the progression cohort is 1.8, 0.5, 0.1, and 2.3 respectively.

Overview of the femoral atlas creation process
The femoral cartilage atlas was generated from the MR images of the right knee of the 30 chosen participants from the incidence cohort using the following steps.

Delineation of the cartilage from the MR knee images
The image on the left top corner in figure 1 shows a slice from the MR images of the knee.Automatic identification of the cartilage from the knee image is a challenging endeavor because of low contrast in some areas and is an area of intense research activity.A robust completely automated algorithm has not yet been validated on a large number of image volumes.However, for this study our focus was to develop the cartilage atlas and a manual segmentation was used.It was critical to be precise in accurately segmenting the cartilage from the entire knee image.For this purpose, two operators (graduate students) in consultation with and verified by an experienced musculo-skeletal radiologist manually segmented the femoral cartilage from the knee images.An iterative approach was taken where the operators consulted the radiologist before and after performing segmentation for small sets of images.Consultations after segmenting each small set ensured high accuracy and reduced fatigue for the operators; this was critical considering the large number of slices for each participant because of the high resolution MR scans (0.3650.3650.70)and 30 datasets for creating the atlas (Tameem et al., 2007(Tameem et al., , 2011)).

Image pre-processing steps
As a first step, a cubic interpolation scheme was applied to the segmented images to achieve an isotropic resolution of 0.3650.3650.365mm 3 for datasets of all 30 participants.A reference subject was selected with age and KL grade close to the average of the cohort as the initial representative of the cohort.All subjects were than aligned to the reference subject.This initial alignment was achieved using a mutual information based affine transform technique available through the FMRIB Software Library (FSL) (Jenkinson & Smith, 2001).The affine transformation corrects for global size and positional differences between a subject and the chosen reference.The affine transformed subject data is then mapped to the reference using an elastic registration technique that is based on the Demons algorithm (Thirion 1998).This mapping results in a 3-dimensional deformation field that spatially maps each voxel in a subject dataset to the coordinate system of the reference.This process is iterative and is computed using the equation 1 (Ardekani & Sinha, 2006).
In equation 1, u n+1 denotes the correction vector field at iteration n+1, G  is a Gaussian filter with variance ,  represents the convolution, scaling factor denoted by C, and T and S denote the reference and transformed image intensities respectively.The displacement u(q) is estimated by the algorithm such that for every voxel location (q) in the reference image (T) it maps the corresponding location in the subject's transformed images (S).To preserve the morphology, the algorithm computes deformation fields through both forward and backward transformations.The Gaussian filter was optimized to a 3x3x3 kernel to make sure the deformation is smooth and registration is more accurate.A mean intensity image is created when the transformed images of the entire group are averaged.Similarly, averaging the deformation fields of all images in the group yields a mean deformation field.Combining this mean deformation field with the mean intensity image creates an image template for the group.In the next iteration the image template produced replaces the reference image and the whole process of affine and elastic transformation for all 30 image set is repeated again.This process continues until there is no substantial change in the deformation field between two consecutive computations.This was observed after 4-5 iterations when creating the cartilage atlas.The final image created is the average shape atlas that converges to the centroid of the population data set (Kochunov et al., 2001(Kochunov et al., , 2005)).The entire process is visualized in figure 1.

Active shape models
Active shape models (ASM) developed in the past (Cootes et al., 1994) have been used extensively to determine the patterns of variability within a group of subjects.Using principal component analysis, ASMs can be created from the deformation fields that are available after the last step in the atlas creation process.We define the deformation field as a 3-dimensional matrix that stores the amount of displacement needed to move a voxel on the atlas to its corresponding voxel location in the subject (this is achieved with the freeform transformation).The analysis was performed in the following steps also discussed in detail in our recent publication (Tameem et al., 2011).The mean deformation value over N subjects is calculated by averaging the deformation field for all subjects over every voxel as shown in equation 2. () Using equation 3 we can calculate the deviation of each subject from the mean value calculated in equation 2.
A covariance matrix of dimension n×n is calculated as show in equation 4.This covariance matrix enables calculating the basis as shown in equation 4. () Finally, the eigenvectors and eigenvalues are calculated by diagonalizing the covariance matrix.Using equation 5 a linear model is created which is our shape model.
In equation 5,  = (1, 2… k) is a matrix of the first eigenvectors, and Ws is a vector of weights called the shape coefficient.The principal component analysis on the deformation fields of all 30 participants results in the lead modes of shape variation.Shape variations are calculated along the two leading modes of variations at 2SD and 3D from the mean image.The shape variance for the femoral cartilage within the group of participants used to create the atlas can be visually represented through the images synthesized using ASMs.

Tensor morphometric analysis to detect local shape changes
Local variations in shape between populations can be estimated using the jacobian values of the deformation fields available from the non-linear warping algorithm.The deformation fields obtained from the last step of atlas creation process provide for a voxel (q) in the atlas reference frame, a three dimensional map of displacement values.For every voxel (q) mapped from the reference image (atlas) to the subject image (i), the displacement values can be broken into three components namely u i , v i , and w i .Each voxel itself can be expressed by its three coordinates x, y, and z.The jacobian J i (q) as shown in equation 6, is defined as the determinant of the gradient mapping function and I is the identity matrix.The jacobian values calculated for each subject in a group are averaged as show in equation 7.
J p denotes to the mean of the jacobian values over all subject in the cohort under consideration.N represents the number of subjects in each population group, and J i (x, y, z) denotes the jacobian value at each voxel location for subject i.When images are warped the calculated jacobian values provide information on the localized volume changes between two populations.Using statistical techniques regions with significant differences between the two populations (the incidence and progression cohort) can be identified on a voxel level.Jacobian maps provide a powerful visual description that highlights these changes.

Statistical analysis
The jacobian values were obtained after aligning each subject in the incidence and progression cohort to the atlas as described in the previous step.Some processing steps such as smoothing the deformation fields using a 4mm full-width-half-maximum Gaussian kernel were performed.Smoothing takes care of any inaccuracies in the registration and ensures a normal distribution.The statistical analysis is performed using the Statistical Parametric Mapping (SPM) application (version 5) in MATLAB (Ashburner, 2009).A twosample t-test voxel based analysis was performed on the dataset between the incidence and progression cohort to determine the local variations in the cartilage morphology between the populations.A false discovery rate (FDR) correction was used to control for multiple comparisons (Benjamini & Hochberg, 1995).FDR is an established method to correct for multiple comparisons in brain voxel-based morphometry.Statistical analysis yields some regions with higher "t" values that indicate areas of larger morphological differences in the femoral cartilage between the incidence and progression group.

Registration and atlas
Accuracy in segmentation and registration are critical components in the atlas creation process.The initial affine transformation correct for positional and size differences.The freeform transformation applied after affine transformation measures true morphological differences.Figure 2 shows the alignment results of one of the slices of a 3D image volume.
The contour shown on the reference image is transferred to the subject image before and after affine and non-linear transformation to confirm accuracy achieved in the registration process.It is fairly evident from the results that the affine transform does a very good job of correcting the positional shape changes, whereas, the non-non-linear transform corrects the local changes.After 3 iterations, no significant change in the deformation field values was observed in successive iterations.The atlas converged to the population centroid after three iterations and the sharp edges seen in the atlas (figure 3) confirm the accuracy of the free form registration.

Visualizing deformation maps
The mean and standard deviation of the femoral cartilage calculated from the magnitude of the 3D deformation fields across 30 subjects are visualized in figure 4 (2D display) and figure 5 (3D display).The regions of large deformations and standard deviation values are seen along the edges (Figure 5).A possible explanation resulting to this finding is the inaccuracies in the manual segmentation process.The trochlea region near the intercondylar notch and in the medial aspect shows higher deformation and standard deviation values (Figure 5).Between the medial and lateral aspects of the cartilage, the medial condyle regions show slightly higher deformation than the lateral regions.

Active shape models from atlas cohort
Synthesized images using principal component analysis based active shape models were created.Using the first 2 dominant eigenvectors the images generated at 2SD and 3SD and the 3D maps are visualized in figure 6.
Fig. 6.Synthesized images along 2SD and 3SD overlaid on the atlas.The first and second row corresponds to images generated along the 1 st and 2 nd eigenmodes respectively.
The synthesized images (along 1 st and 2 nd eigenvectors and different standard deviations) are overlaid on the atlas.Regions where the atlas and the synthesized images completely overlap are shown in yellow.Regions where either the atlas or the synthesized images extended beyond each other are shown in red and green respectively.It should be noted that the shape of the atlas remains the same and the extension of synthesized images beyond the atlas changes for different standard deviations.Looking closely at the synthesized shape images at the positive and the corresponding negative SD it is observed that the areas visualized in red (and green) in one image are visualized as green (and red) in the corresponding image.Also, the magnitude of shape difference is evidently higher in at 3SD than at 2SD|.

Statistical analysis
Statistical analysis performed in SPM results in 't' values that are obtained after comparing the jacobian values for incidence and progression cohorts.The regions with significant differences are overlaid on the atlas and can be clearly visualized in the 3D map shown in figure 7. Regions such as the medial weight-bearing region and the lateral posterior condyle exhibit significant differences between the incidence and progression cohorts.

Discussion
Extending the previous work of Cohen et al. and Carballido-Gamio et al. (Cohen et al., 2003;Carballido-Gamio et al., 2008), the focus of this study was create a femoral cartilage atlas from high-resolution 3D MR images made available by the OAI through version 0.A.1.The atlas-based approach to study localized morphological changes has been well established for brain and we extend this to study localized changes in cartilage morphology in osteoarthritis.Figure 2 confirms the accuracy of the registration steps involved in the atlas creation process (affine and freeform) by overlaying the contour traced on the atlas over the subject image obtained before and after the two registration steps.In addition to this visual confirmation, 3D residual distance maps have confirmed that the nonlinear algorithm used here provides registration accuracy within 1-2 pixels.
Figure 3 shows the 2D slices and 3D volume of the femoral cartilage atlas that was created from 30 Caucasian male subjects from the incidence group.The accuracy of the registration techniques (both affine and freeform) can be confirmed by the sharp edges observed in the 3D image.Several non-linear deformation algorithms are available to warp an image volume to a target volume.The accuracy of the algorithm is usually verified by the accuracy of the registration.However, this still does not guarantee that the jacobian values calculated from the deformation fields reflect the actual local volume changes.Rohlfing et al. has confirmed, using synthetic image volumes and deformations, that the demons algorithm used in the current work does provide accurate estimates of local volume changes (Rohlfing, 2006).The deformation values obtained for the 30 subjects used to create the atlas are averaged and the voxel-wise mean and standard deviation values are visualized in figure 4 and 5.It is observed that the overall deformation values are quite small (in the order of ~0.4 mm) and hence exhibit very less variation within the population used to create the atlas.However, higher values are observed around the edges and that can be attributed to registration errors.The low variability seen in the current atlas can be credited to the fact that the participants selected to create the atlas had to meet strict selection criteria of age, gender, race and low K&L score.Past studies found the inter-subjects variability across morphological parameters such as cartilage thickness and volume to be significantly higher (Eckstein et al., 2001).In our study the inter-subject variability was minimized by using the affine transformation as the first step in the atlas creation process and accounts for global size changes.In figure 5, higher mean and standard deviation of the deformation values are realized on the patellofemoral compartment of the trochlear (femoral) side.Previous distal femoral bone studies have indicated that shape changes in the intercondylar notch when comparing normal against diseased subjects as a discriminant of OA (Shepstone et al., 2001) and hence has there is a possibility of a potential link between the findings of this study to what has been observed in the past.
The active shape models generated using principal component analysis along first two eigenmodes and ±2SD and ±3SD are shown in figure 6.These models display small morphological variability that is in agreement with the findings using the deformation magnitude values.Integrating active shape models with segmentation algorithms can provide a-prior information about cartilage shape variability (Duchesne et al., 2002).Also, active shape indices derived can be used to discriminate between cohorts or classify single subjects.
Results from the voxel based statistical analysis performed using SPM in Matlab as shown in in figure 7 display specific weight bearing regions such as the lateral side of the intercondylar notch and medial posterior femoral regions displaying significant morphological differences between the incidence group and the progression group.These findings are consistent with other studies.E c k s t e i n e t a l .s h o w e d t h a t c a r t i l a g e degradation is significantly higher in knees with radiographic evidence of OA in comparison to healthy knees and knees with no radiographic evidence of OA.Subjects were grouped as non-exposed controls (n = 112), calculated K&L grade 2 (n = 310), calculated K&L grade 3 (n = 300), and calculated K&L grade 4 (n = 109).Regional and subregional thickness values calculated from coronal MR images showed that there was significant difference seen in the weight bearing sub-regions of the femorotibial cartilage in subjects with K&L grade of 4 (Eckstein et al., 2010).Bredbenner et al. recently investigated the potential of using statistical shape modeling as a tool to detect onset of OA.The authors use SSM to effectively characterize the variability in the subchondral bone region in femur and tibia.They also show the potential of combining SSM with rigid body transformation to distinguish subjects at risk and no risk of developing OA (Bredbenner et al., 2010).
As an extension to our current work, we have created statistical atlases of the femur and tibia bone.We aim to use these atlases to characterize morphological variations in the bone in a well-defined group and use it as a predictor for osteoarthritis and osteoporosis.For this study MR images of the knee were acquired on 10 normal subjects and the femur and tibia were segmented from the sagittal MR images.The atlas of the femur and tibia are created as outlined in section 2.2.The accuracy of registration is shown in figure 10 where the contour traced on the atlas is overlaid on the test subject and the remaining images.
The atlas for femur and tibia are visualized in figure 9.The sharp edges and the thin growth plate seen in the cross-sectional images is a visual confirmation of extremely accurate alignment during the atlas creation process.We present some initial results and in future hope to utilize the rich data set available from OAI and conduct a detail study to investigate morphological variation in the bone using the atlas based approach.

Conclusion
This chapter is focused on the atlas-based approach for automated image analysis for cartilage and bone.The primary area of application is to characterize and monitor progression of osteoarthritis.In this chapter, the atlas approach is applied to the study of morphological differences between population cohorts at different stages of osteoarthritis.
The second area of application is a preliminary study of femoral and tibial bone atlas with the aim of identifying morphological differences in bone in matched cohorts with osteoarthritis.The atlas approach is not limited to identification of morphological differences (between cohorts) or changes (in longitudinal studies).It is readily extendable to the analysis of parametric image datasets such as T2, T1rho and dGEMRIC.

Fig. 1 .
Fig. 1.An illustration of various steps involved in the atlas creation process: starting from the initial data selection, manual segmentation, affine registration, and freeform registration.Iterating the freeform registration and updating the template results in convergence to the centroid of the cohort.

Fig. 2 .
Fig. 2. Image showing registration accuracy achieved during the atlas creation process.(From left to right): The contour traced on the reference image is overlaid on the subject image, subject image after affine transform, and subject image after freeform transformation.

Fig. 3 .
Fig. 3. Figure illustrates the 2D slices of the atlas image in the left column and a 3D image in the right column.The sharp edges displayed on the atlas image are evidence of high registration precision (S-superior, L-lateral, and M-medial).

Fig. 4 .
Fig. 4. Two-dimensional maps of the mean and standard deviation values of the magnitude of the 3D deformation vector fields of the 30 subjects used to create the atlas.The values are displayed over several cross-sectional slices in the entire atlas volume.

Fig. 5 .
Fig. 5. Three-dimensional map of the mean and standard deviation of the magnitude of the 3D deformation vector fields across the 30 subjects used to create the atlas.

Fig. 7 .
Fig. 7. Results from the statistical analysis overlaid on the atlas.Regions with most difference between the incidence and progression cohort is shown on the 3D map.The regions in white indicate no change.

Fig. 8 .
Fig. 8. Selection on the reference image is overlaid on the subject before and after affine and freeform registration.For Row 1: Femur and Row 2: Tibia

Table 1 .
Table 1 provides the detailed demographics of all participants.Demographics of participants from incidence and progression sub-cohort The jacobian values that are calculated are always positive.A value of 1 indicates no volume change, values greater than 1 represent positive volume change, and values lower than 1 represent a negative volume change.