Registration of Magnetic Resonance Image Series for Knee Articular Cartilage Analysis: Data from the Osteoarthritis Initiative.

OBJECTIVE: Although conventional radiography is used to assess osteoarthritis in a clinical setting, it has limitations, including an inability to stage early cartilage degeneration. There is a growing interest in using quantitative magnetic resonance imaging to identify degenerative changes in articular cartilage, including the large multicentered study, the Osteoarthritis Initiative (OAI). There is a demand for suitable image registration and segmentation software to complete this analysis. The objective of this study was to develop and validate the open source software, ImageK, that registers 3 T MRI T2 mapping and double echo steady state (DESS) knee MRI sequences acquired in the OAI protocol.
METHODS: A C++ library, the insight toolkit, was used to develop open source software to register DESS and T2 mapping image MRI sequences using Mattes's Multimodality Mutual information metric.
RESULTS: Registration was assessed using three separate methods. A checkerboard layout demonstrated acceptable visual alignment. Fiducial markers placed in cadaveric knees measured a registration error of 0.85 voxels. Measuring the local variation in Mattes's Mutual Information metric in the local area of the registered solution showed precision within 1 pixel. In this group, the registered solution required a transform of 56 voxels in translation and 1 degree of rotation.
CONCLUSION: The software we have developed, ImageK, provides free, open source image analysis software that registers DESS and T2 mapping sequences of knee articular cartilage within 1 voxel accuracy. This image registration software facilitates quantitative MRI analyses of knee articular cartilage.

Introduction

Plain radiography is the traditional imaging standard used to diagnose osteoarthritis (OA). Disease progression can be tracked following serial radiographs, assessing changes in joint space width, and the appearance of other OA hallmarks including osteophytes.[1] The joint space narrowing observed on radiographs does not occur until there has been significant cartilage damage, typical of the late-stage disease.[2] Early OA is thought to involve microstructural damage of cartilage without significant volume loss. The inability of radiography to directly image cartilage and assess the early stages of cartilage damage limits its usefulness as a diagnostic and disease-staging tool.

There is a growing interest in using quantitative magnetic resonance imaging (qMRI) to identify the degenerative changes of OA in cartilage. qMRI provides a noninvasive means to assess articular cartilage morphology and molecular biochemistry in three dimensions. Attempts to monitor OA with MRI have focused on documenting morphometric changes including volume loss in OA[3-5] and molecular changes including anisostropy of collagen fibers and water and proteoglycan content.[6-9] One large multicentered study, the Osteoarthritis Initiative (OAI), has collected qMRI cartilage data using double echo steady state (DESS) and T2 mapping sequences. DESS imaging provides volumetric information, whereas quantitative T2 mapping probes aspects of water content and collagen fiber orientation.[7-10] To access the information contained in the T2 mapping data, there is a need for freely available software to register the T2 maps to higher resolution morphologic sequences, such as DESS, on which accurate cartilage segmentation can be performed.

To complete quantitative analysis on T2 maps, individual values of each pixel of cartilage needs to be extracted from the image. This can be accomplished by cartilage segmentation. In large patient cohorts, such as the OAI, hand segmentation of the cartilage becomes unreasonable because of the large time commitment required. Automated cartilage segmentation for T2 maps is challenging because of the poor signal contrast between cartilage and soft tissue. Many groups have previously developed automated segmentation software for DESS images because of the increased cartilage soft tissue contrast.[11-13] These algorithms have segmented the entire cartilage volume and have not differentiated different layers of cartilage. The automated segmentation software for DESS images could be applied to T2 maps to extract the quantitative T2 values if the DESS and T2 map sequences were registered.

In a typical registration problem, collected pixel intensities are derived from the same type of signal. Consequently, metrics in the registration algorithm used to compare the degree of alignment between the image sets can be derived from direct comparison of pixel intensities in corresponding locations of each image. In multimodality registration, such as between DESS and T2 maps, pixel intensities cannot be directly compared because the same anatomic structure on DESS and T2 sequences do not necessarily possess the same pixel intensities. Therefore, a more sophisticated metric is required for multimodality registration, typically a mutual information framework.[14,15] In this work, Mattes’s Mutual Information metric was implemented for multimodality registration.[16]

Here, we develop and describe open source software that allows rigid registration of DESS and T2 mapping sequences. Furthermore, we quantified the accuracy and precision of the registration process. This software can be used with manual, semiautomated, or automated segmentation of DESS images to accurately segment articular cartilage on T2 mappings. The software is offered under the BSD open source software license and is available with example data sets included as supplemental data for free download (www.imagek.org or at the permanent archived site: http://hdl.handle.net/10380/3355).

Methods

MRI Acquisition

Three-dimensional DESS with water excitation images and T2 mapping images were acquired using sequences approved for the National Institutes of Health–sponsored Osteoarthritis Initiative study.[17] MRI of the knee joint was performed on a 3.0-T Siemens whole body MAGNETOM Trio 3T scanner (Siemens, Erlangen, Germany) using a standard extremity coil. For high-spatial-resolution 3D DESS imaging before registration,[10] a total of 160 sections were acquired with a field of view (FOV) of 14 cm (matrix 384 × 307) with an in-plane spatial resolution of 0.365 × 0.456 mm and a slice thickness of 0.7 mm with an acquisition time of 11 minutes. For sagittal 2D dual-echo fast spin echo (FSE) sequence for mapping T2 relaxation time, TR was 2700 ms and 7 echo images with TE ranging from 10 to 80 ms were acquired with matrix of 384 × 384, in-plane resolution of 0.313 × 0.313 mm, FOV of 12 cm, acquisition time 11 minutes, and slice thickness of 3 mm.

Registration

Multimodality rigid registration of DESS and T2 images were completed using Mattes’s Mutual Information metric across a standard framework (). The initial description is applied using a deformable model.[16] Here, we assumed the bodies were rigid. The registration software was built using the insight toolkit, a C++ open source, cross-platform image analysis library (www.itk.org).[18] DESS sequences were rotated and translated through three dimensions in the registration process (the “moving” image set) because of their higher and more isotropic resolution as compared with the fixed T2 mapping sequences (the “fixed” image set). Registration between the T2 mapping multiple echo sequences in the multi-echo spin echo (MESE) set was unnecessary as the echo images are not acquired serially but as a single interleaved acquisition. Only the first echo of the MESE set was used for registration since each echo sequence was assumed to be registered to the other echo sequences in the MESE set. In the DICOM header file, the DICOM tag (0018, 0081) was used to extract the first echo image series from the entire MESE set. The tibia and femur were registered simultaneously. Images were transformed through physical space using a three-dimensional rigid verser transform where rotation is represented using a unit quaternion and translated by a vector. During the transform, linear interpolation was used to estimate pixel values during resampling when an exact pixel location on the image grid was not available from the sampled point in physical space (i.e., a nongrid position). The registration search space is large, across 6 degrees of freedom. As a result, a specialized gradient decent optimizer is used to define the criterion for the transform parameters to move efficiently through the registration search space. After the transform, the mutual information metric is used to assess the degree of alignment between the two images, and the process is repeated until a maximum degree of overlap has been achieved. Parameters were optimized to ensure robustness of accuracy. Registration was completed on an HP590t series computer with a 3.33 MHz Intel core i7-980X six core processor, 24 GB SDRAM, 2TB 7200 rpm SATA hard drive, and a Windows 7 64-bit operating system (Hewlett Packard, Pala Alto, CA). After registration, the DESS images were interpolated to contain the same image resolution as the initial T2 mapped images.

Figure 1.

A schematic outline of the registration framework is shown.[13]

Registration Validation

Registration error was determined by comparing the center of mass of multiple MR lucent fiducial markers on corresponding T2 mapping and DESS sequences. In three cadaveric knees (LifeLegacy Foundation), MR lucent beads were fashioned from 4% agar doped with 2 mM Gd-DTPA and placed in separate drill holes in the femoral condyle and tibial plateau of each knee for a total of four to five markers in each knee (Magnevist, Berlex Imaging, Wayne, NJ). Markers had an approximate diameter of 1 cm. Each knee was imaged as described above and registered. Identical MR sequencing protocols were used as described for the OAI.[17] No attempts to significantly alter the motion of the leg were made between sequences. The coordinates for the centers of each marker were calculated in the DESS images and compared with the corresponding center of mass for the marker on the T2 mapping images by hand segmenting the bead in both images and binary eroding the segment until one pixel remained. Density distributions were estimated using two-dimensional Gaussian kernels.

Registration precision was determined by comparing the Mattes’s Mutual Information metric value at the registration solution to a series of given image displacements to introduce known registration error. The second release of images for the OAI was used. Data used in the preparation of this article were obtained from the OAI database, which is available for public access at http://www.oai.ucsf.edu/. Specific data sets used from the OAI were kXR SQ reading (BU), version 0.5, for measurement of the centrally read Kellgren and Lawrence (KL) grade scores (variable identification: V00XRKL). Three populations were extracted from the OAI: the entire nonexposed control cohort (n = 122), the standard subset of the progression cohort (n = 160) released by the OAI, and all subjects in the progression cohort with a KL grade of 4 (n = 142). All DESS and T2 mapping sequences from these groups at the initial baseline time point were registered as described above. The transform for the registered solution between the DESS and T2 mapping sequences in each patient was defined as the coordinate origin (0, 0, 0). The DESS sequence was then translated a given distance across the T2 mapping sequence to create a known degree of misalignment. The Mattes’s Mutual Information metric value was recalculated at the new translated coordinate and compared with the metric value at the initial coordinate origin. This was repeated across a range of translations extending across 5 voxels in each direction along all three axes. During registration, it is possible to rotate the image across three degrees of freedom. As such, at each translated coordinate, the DESS sequence was rotated across all three axis in all possible rotation combinations in the range (−1/90·π, 1/90·π) or 2° in each direction in increments of 1/720·π (0.25°). The rotation with the largest Mattes’s Mutual Information metric value (indicating the lowest degree of alignment) at the given translation was selected as the representative value for that transformed coordinate for visualization purposes. Isometric contour maps of the normalized percent variation in the Mattes’s Mutual Information metric value as a function of translation were generated using the statistical software package R.[19]

Statistics

Data are expressed as mean ± confidence interval, except where noted. Direct comparisons between two populations were made using an unpaired, two-tailed Student’s t-test. Statistical significance was determined if P < 0.05. Multiple-group comparisons were made using two-way ANOVA, using the Student–Newman–Keuls pairwise comparison to determine significance levels. All statistical tests and analysis were completed using the statistical software package R (www.r-project.org).[19]

Results

Registration Error

The accuracy of registration was validated using two separate approaches. Accuracy of registration was quantitatively measured using MR opaque fiducial markers. DESS and T2 mapping sequences were collected on three cadaveric knees with MR lucent fiducial markers placed in the femoral condyle and tibia plateau. Registration was performed, and the center of mass of each marker was compared with their retrospective location on DESS and T2 mapping sequences. The registration error across the x- and y-axis was measured at 0.85 ± 0.3 voxels (SEM; ). In the z-axis, alignment was identical; however, given the large difference in resolution of the z-axis (0.7 mm for DESS vs. 3 mm for T2), the precise degree of pixel alignment could not be determined below the z-axis resolution of the T2 image sets (3 mm). The addition of the fiducial markers to the image sequence could potentially improve registration accuracy by providing additional cues for the registration process. To verify that the addition of the fiducial markers was not altering accuracy of the registration methods, registration was completed with and without the image regions that contained the fiducial markers. Comparing registration transforms between these two groups resulted in a difference in total translation of less than 0.27 ± 0.1 voxels and rotation was less than 0.03 ± 0.01 degrees of rotation of the transformed images and was not a statistically significant difference when compared with the inclusion of fiducial markers.

Figure 2.

Registration of T2 mapping and double echo steady state (DESS) sequences using multimodality registration is accurate. Qualitative accuracy is shown using a checkerboard layout where T2 and DESS sequences are shown overlaid in alternating squares producing a checkerboard appearance. An unregistered (A) and the registered image (B) are shown. For comparison purposes, the anterior aspect of the femoral condyle of each sequence in each image is outlined. (C and D) Representative images of registered T2 and DESS sequences show the precision of alignment. The arrow notes the posterior cruciate ligament. Scale bar in bottom right of (D) represents 20 voxels. (E) The femoral and tibial cartilage interface outlined in (D) is enlarged to show alignment of fine structures after registration. Scale bar represents 8 voxels. Multiple fiducial magnetic resonance opaque markers were placed in three separate cadaveric knees and the center of mass of each marker between both the T2 and DESS sequences were compared to quantitate accuracy. (F) The area of disagreement between registered and ground truth pixel locations of both the x- and y-axis are represented as a Gaussian kernel density estimate on the two-dimensional histogram. The mean value was 0.85 ± 0.3 voxels (SEM).

Second, a checkerboard layout provided a tool to visually assess the accuracy of registration. As compared with the initial unregistered overlapped images (), postregistration demonstrates subpixel accuracy (). A series of image slices across the same individual knee show subpixel alignment across the femoral and tibial articular cartilage and ligamentous structures including the posterior cruciate ligament ( and ). Examination of the femoral and tibial cartilage interface shows subpixel agreement (). This accuracy is present throughout the entire three-dimensional sequence (Suppl. Fig. S1). The provided software includes an overlapped checkerboard layout of the image of the two registered images to quickly visually assess registration accuracy on any image set.

Registration Precision

To verify that the severity of OA did not alter registration accuracy or precision, a series of image groups from the OAI were compared. The DESS and T2 mapping sequences of patients from the entire nonexposed control cohort (n = 122), the standard subset of the progression cohort (n = 160) released by the OAI, and all subjects in the progression cohort with a KL grade of 4 (n = 142) at the initial baseline time point were registered. The total translation and rotation required to register the DESS images to the T2 mapped images had a mean translation of 56.3 ± 2.8 voxels and rotation of 1.3 ± 0.5°. In the nonexposed control group, the mean translation and rotation was 57.9 ± 4.0 and 0.49 ± 0.1 voxels, respectively, with an average KL of 0.1 ± 0.1. In the progression cohort, the mean translation and rotation was 59.9 ± 5.7 and 1.7 ± 1.1 voxels, respectively, with an average KL of 2.0 ± 0.2. In the cohort of subjects with a KL grade of 4, the mean translation and rotation was 50.9 ± 3.7 and 1.5 ± 0.7 voxels, respectively. Only the rotation in the nonexposed control group was statistically significant when compared with the other groups. All other rotations and translations between these three groups were not statistically significant.

Registration precision and accuracy was measured using displacement by comparing the Mattes’s Mutual Information metric value at the registration transform solution to other possible transforms in the local neighborhood in each of these three groups. The Mattes’s Mutual Information metric is a mutual information similarity criterion measuring the degree of alignment between two overlapping images. The registration solution is determined by the minimal Mattes’s Mutual Information metric by transforming the DESS image across six degrees of freedom (three degrees of translation and three degrees of rotation) as compared with the stationary T2 mapping sequence. For each subject, the registration solution was calculated between the DESS and T2 mapping sequence. A series of given displacements in the local neighborhood of the registration solution was then introduced by translating and rotating the image. This introduced a given, known registration error from the determined registration solution. The Mattes’s Mutual Information metric between this new transform and the registration solution was compared to identify the precision of the method and to determine if an accurate minimum value had been identified. Registration of each subject was completed within 10.1 ± 0.5 minutes.

The minimal value existed in an approximate 1 voxel neighborhood around the coordinate origin. Outside of this area, sufficient variation in the Mattes’s Mutual Information metric exists at different translations to differentiate between the determined registration solution (the minimum Mattes’s Mutual Information metric value) and less accurate transforms (-). The 0% contour line had a range of 0.98 ± 0.4, 0.78 ± 0.5, and 0.71 ± 0.5 voxels for the nonexposed control, progression, and for the KL grade of 4 group, respectively. The z-axis is not shown given the large difference in resolution between DESS and T2 imaging modalities along this axis as discussed above.

Figure 3.

Registration of T2 maps and double echo steady state (DESS) sequences are precise. The multimodality metric used in registration was measured across a series of translations and rotations in the neighborhood of the calculated optimum solution. Three groups were tested: (A) the Osteoarthritis Initiative (OAI) nonexposed control group, (B) the standard release subset of the OAI progression cohort, (C) and the entire population from the progression cohort with a KL grade of 4. For each group, images were registered, transformed across a local range of points, and the percent difference between the Mattes’s Mutual Information metric value was calculated by comparing the Mattes’s Mutual Information metric at each transformed point to the computer determined registration transform at the origin (0, 0). Incremental isolines indicate an increasing error in registration. The isoline of 0% has an approximate 1 voxel radius from the registration solution. Each contour line represents a 2.0% change in the Mattes’s Mutual Information metric value. Each labeled point noted on the contour map indicates a separate translation that was tested in 0.25 voxel increments in each axis across all possible combinations.

Discussion

To complete quantitative analysis of T2 mapped images, segmentation of the knee cartilage on each T2 mapped sequence is required. In large patient cohorts, such as the OAI, this is a challenging problem given the time commitment necessary to complete hand segmentation. Various groups have developed automated or semiautomated cartilage segmentation software for DESS images. T2 mapped images have low contrast of the cartilage–soft tissue boundary, making automated segmentation extremely challenging. DESS sequences have an increased signal contrast between cartilage–bone and cartilage–soft tissue interfaces. As a result, there are multiple approaches reported in the literature to segment knee articular cartilage on DESS MRI images.[11,12] This software could be applied to T2 maps if the two image sequences could be registered so that the segmentation created for the DESS image could be applied to the T2 map. It should be noted that as a rigid registration method was used this software has not been validated for application to subjects over long time periods where anatomic deformation of the bone, soft tissues, and cartilage may have occurred.

Here, we describe a strategy to align MRI T2 maps and DESS sequences collected from the same knee using multimodality rigid registration. Multimodality registration allows the segmentation of DESS sequences to be applied to T2 maps. After a DESS sequence has been segmented, the cartilage mask can be applied to the corresponding T2 map sequence. At the bone–cartilage interface segmented on DESS images, the non-fat-suppressed MESE sequences have a low signal when compared with fat-suppressed MESE sequences, referred to as the signal void zone. Multiple postprocessing options exist in dealing with extracting quantitative T2 maps using DESS segmentation that contain the signal void zone, including using a threshold between a minimum and maximum physiologic T2 values as a filter. The signal void zone at the subchondral bone interface in these sequences has a thickness of approximately 2 pixels and can be ignored if it does not significantly alter T2 value calculations.

Three other registration methods have been applied to register knee articular cartilage. These methods were based on registration of a single MRI sequence across multiple time points using deformable and rigid registration.[20-22] Here, we describe and validate a method for multimodality rigid registration. Registration error in these previous methods was reported at approximately 1 voxel[20] and at 25% of a voxel volume.[22] On a registration problem of increased difficulty from the perspective of using multiple MR modalities, where these previously reported methods cannot be applied, we report an error of less than 1 voxel.

The OAI cohort was used to generate a large population of sequences to define the precision and confirm the accuracy of this method. By definition of the principles of mutual information, homogenous regions of the DESS image set directly map to homogenous regions of the T2 map image set.[14,15,23,24] Consequently, by definition, the optimum Mattes’s Mutual Information metric (i.e., the minimal value) predicts the best possible registration transform between the two image sets. This does not necessarily guarantee perfect alignment between the two sequences but does represent the highest possible accuracy that can be theoretically obtained. The optimal registration transform was identified within a 1 voxel neighborhood of the calculated solution, supporting our conclusion of registration accuracy. Furthermore, the variation in Mattes’s Mutual Information metric values in the local neighborhood beyond 1 voxel of the optimum solution confirmed the method’s precision.

Misalignment and error in the registered images is a result of subtle movements of the subject between acquisitions, changes in field of view, and differences in prescribed slab geometry. As a result, the images cannot be registered by using the DICOM header files because they do not contain the images geometric information to align the two sequences. Direct overlap of the two sequences prior to registration results in obvious, unacceptable error (). Average misalignment between sequences was measured at 56 voxels and 1 degree of rotation. The OAI acquisition protocol is stringent; MR technicians are educated during acquisition to only translate and rotate the field of view between acquiring the DESS and MESE sequences. This likely explains the low degree of rotation needed to register sequences. Any misregistration is not a function of the time between sequencing or of changes in anatomic geometry, and a deformable registration model is not necessary. A limitation in the technique includes simultaneous registration of the femur and tibia as a rigid structure, instead of registering the tibia and femur separately. These two articular surfaces may develop slight differences in their relative positions after possible motion through the joint in the brief time period between DESS and MESE sequences.

The method is robust. The large size of the OAI cohort provides a large series of repeated measurements, at different institutions, with different operators, and operational repositioning between measurements that demonstrates the robustness of the software. The process is completely automated, allowing rapid and easy registration of large data sets such as the OAI.

The approach to registration described in this article can be used to assist segmentation of articular cartilage in T2 maps from DESS images with any segmentation method. This strategy may also be applied to register DESS or fat-suppressed T1-weighted spoiled gradient echo images with T1ρ or dGEMRIC maps. The accuracy of our software has not been validated with these sequences; however, the included “checkerboard layout” software allows fast and rapid validation of such application.

The open source software movement in image analysis is based on reproducible results. Complete access to the code and executables with similar data files used by the authors is freely available to the community to test the reproducibility of the findings in this article and to maximize the potential of public image databases, such as the OAI. Overall, we have demonstrated that the two image sequences, DESS and T2 maps, can be registered within 1 voxel accuracy and precision.