- Research article
- Open Access
Development of a morphology-based modeling technique for tracking solid-body displacements: examining the reliability of a potential MRI-only approach for joint kinematics assessment
BMC Medical Imaging volume 16, Article number: 38 (2016)
Single or biplanar video radiography and Roentgen stereophotogrammetry (RSA) techniques used for the assessment of in-vivo joint kinematics involves application of ionizing radiation, which is a limitation for clinical research involving human subjects. To overcome this limitation, our long-term goal is to develop a magnetic resonance imaging (MRI)-only, three dimensional (3-D) modeling technique that permits dynamic imaging of joint motion in humans. Here, we present our initial findings, as well as reliability data, for an MRI-only protocol and modeling technique.
We developed a morphology-based motion-analysis technique that uses MRI of custom-built solid-body objects to animate and quantify experimental displacements between them. The technique involved four major steps. First, the imaging volume was calibrated using a custom-built grid. Second, 3-D models were segmented from axial scans of two custom-built solid-body cubes. Third, these cubes were positioned at pre-determined relative displacements (translation and rotation) in the magnetic resonance coil and scanned with a T1 and a fast contrast-enhanced pulse sequences. The digital imaging and communications in medicine (DICOM) images were then processed for animation. The fourth step involved importing these processed images into an animation software, where they were displayed as background scenes. In the same step, 3-D models of the cubes were imported into the animation software, where the user manipulated the models to match their outlines in the scene (rotoscoping) and registered the models into an anatomical joint system. Measurements of displacements obtained from two different rotoscoping sessions were tested for reliability using coefficient of variations (CV), intraclass correlation coefficients (ICC), Bland-Altman plots, and Limits of Agreement analyses.
Between-session reliability was high for both the T1 and the contrast-enhanced sequences. Specifically, the average CVs for translation were 4.31 % and 5.26 % for the two pulse sequences, respectively, while the ICCs were 0.99 for both. For rotation measures, the CVs were 3.19 % and 2.44 % for the two pulse sequences with the ICCs being 0.98 and 0.97, respectively. A novel biplanar imaging approach also yielded high reliability with mean CVs of 2.66 % and 3.39 % for translation in the x- and z-planes, respectively, and ICCs of 0.97 in both planes.
This work provides basic proof-of-concept for a reliable marker-less non-ionizing-radiation-based quasi-dynamic motion quantification technique that can potentially be developed into a tool for real-time joint kinematics analysis.
Visualization of skeletal elements is central to three-dimensional (3-D) kinematic analysis of joint motion. Indirect methods based on tracking surface landmarks (using reflective markers attached to the skin surface) within a calibrated volume (stereophotogrammetry) can contain artifacts (errors of transformation) due to integumentary displacements relative to actual skeletal motion [1–5]. Direct visualization of bony elements during joint motion are typically accomplished via fluoroscopy or cineradiography. Unfortunately, both of these techniques require the use of ionizing radiation, and outcomes from these techniques are restricted mostly to two-dimensional (2-D) analyses as the majority of these systems use single-plane imaging [6, 7]. Emergence of the roentgen stereophotogrammetry (RSA) technique has enabled in-vivo measurement of complex 3-D skeletal kinematics from a series of radiographs acquired with biplanar, orthogonal fluoroscopy [1, 8, 9]. Although this technique is accurate, it commonly requires surgical implantation of markers in bones [2, 9–12], although model-based RSA techniques have recently begun to appear in the literature [13–15].
Recording a series of joint-motion images using x-ray fluoroscopy and then manually superimposing 3-D models of the same skeletal elements to match corresponding outlines in the x-ray images has been used to quantify in-vivo joint motion. [7, 16–18]. More recently, Gatesy et al. reported using the scientific rotoscoping (SR) motion analysis technique, which involves biplanar fluoroscopy to image skeletal movements, creation of 3-D models of joint skeleton from high-resolution computed tomography (CT) scans, followed by model-to-image matching and registration (rotoscoping) performed over several frames of images yielding skeletal motion animation and 3-D kinematic data [6, 19–21]. SR was developed from the X-Ray Reconstruction of Moving Morphology (XROMM) motion quantification technique, which tracks implanted markers digitized in biplanar fluoroscopic images captured within a calibrated imaging volume, instead of utilizing the model superimposition technique . Though accurate, both SR and XROMM techniques require corrections of geometrical distortions in images prior to the animation [11, 12, 22–24]. While x-ray-based motion analysis techniques like SR, XROMM, and RSA are clearly novel and advanced, their translation to clinical research (i.e., human subjects research) has been limited due to health-related concerns associated with the radiation exposure [25–29].
From a clinical research perspective, understanding in vivo skeletal motion is of interest to both scientists and clinicians [30–32]. More specifically, x-ray-based diagnostic imaging techniques measuring human inter-vertebral displacements have focused mostly on imaging the spine at static end-of-range positions [32–40]. However, qualitative and quantitative assessments of spinal motion have been enhanced by quantitative radiographic techniques that track displacements of pre-assigned coordinate points of specific anatomic locations on orthographic spinal images and by real-time joint-motion evaluation with XROMM-like techniques (using CT/magnetic-resonance-imaging-based 3-D models) and RSA (with per-operative implanted vertebral markers) in human subjects [22, 28, 38, 41–47]. Regrettably, these approaches still require exposure to ionizing radiation and, at times, require marker implantation on the bones.
Magnetic resonance imaging (MRI), when used for quantifying inter-vertebral motion, has mostly been restricted to the analysis of end-of-range sagittal-plane displacements [48–50]. However, dynamic cine-phase contrast (cine-PC) or fast-phase contrast (fast-PC) imaging with ultra-fast gradient echo sequences have been employed for evaluating joint kinematics, especially in ankle, knee, or shoulder motion [51–57]. The main approach for these techniques has been the use of pulse sequences that permit volume extraction from full 3-D motion datasets at selected time points along the range of motion (ROM). However, these techniques can be time-consuming. Additionally, the use of cine-PC sequences require a repeated, cyclic, velocity-controlled motion to be performed at the joint of interest during scanning to make the motion synchronized with velocity-encoded motion capture [57, 58]. Also, these images have low resolution and may present motion artifacts [56, 58, 59]. More recently, the combined use of segmented 3-D anatomical models (from high resolution, ~15 mins duration, static axial scans) registered to low resolution, volumetric images acquired at different joint positions using high speed (~40 sec) T1 sequences has been reported . Although such techniques acquire multi-position data with much greater speed, the segmentation of these low-resolution images still require multi-slice images of the experimental quasi-dynamic joint positions. Accordingly, recent advancements in these methods have focused on the acquisition of faster and fewer slices of joint motion (without compromising image resolution) for model-to-image registration and without reducing the accuracy of the technique (time-accuracy tradeoff).
Currently, no modeling techniques exist for quantification of inter-vertebral joint displacements using single-plane or orthogonal magnetic-resonance (MR) image templates for 3-D model registration. Accordingly, our long-term goal is to develop a 3-D model-based technique that permits fast dynamic MR imaging of the human lumbar spine using an open-bore weight-bearing musculoskeletal MRI. Our study focuses on the lumbar spine as low back pain (LBP) is one of the most common reasons for seeking medical care world-wide and accounts for over 3.7 million physician visits per year in the United States alone [60–64]. As such, LBP is arguably one of the most debilitating and costly health disorders, and the development of technologies to aid scientists and clinicians in better understanding the etiology of LBP—as well as in monitoring the effects of therapeutic interventions— is desperately needed.
As a first step towards our long-term goal, we present in this article our initial research and development findings for an MRI-only protocol involving imaging (using a standard T1 and a fast contrast-enhanced MRI sequences), a series of pre-determined displacements between solid-body models, 3-D models (segmenting T1 weighted axial scans), and a morphology-based rotoscoping strategy for animation and quantification of the displacements. The use of the contrast-enhanced sequence will allow us, firstly, to test the feasibility and reliability of its use as a fast imaging tool and secondly, to compare its outcome with that of the standard high-resolution T1 images. The feasibility and reliability of this MRI-based technique is discussed here, and we anticipate further developing this technique into a motion-assessment tool for the lumbar spine and other di-arthrodial joints.
General overview of the experimental design
The experiment involved scanning a pair of wooden cubes placed at pre-determined positions (displacement trials) relative to each other in a custom-calibrated coil of an open-MRI system (0.3 Tesla; Esaote G-scan Brio, Genoa, Italy). Additional axial images of the solid cubes were acquired and segmented using the AVIZO software (Hillsboro, OR, USA) to create 3-D virtual models of the cubes. Next, the MR images of the displacement trials and the 3-D cube models were transferred into an animation software (AutoDesk MAYA, San Rafael, CA, USA); and animations of these displacement trials were performed to quantify the relative motion incurred by the solid bodies. The technique involved four principal steps (Fig. 1a). First, the imaging volume of the MRI coil was calibrated using a custom-built grid (Fig. 2a). Second, 3-D models were segmented out from axial scans of the solid-body cubes (Fig. 2b ii-v). Third, the solid bodies were positioned at pre-determined displacements relative to each other in the MRI coil and scanned (Fig. 2b i); and the digital imaging and communications in medicine (DICOM) images were pre-processed into gray-scale TIF format. Fourth, these images were imported into the animation software using calibration data acquired from the grid used in the first step. These images were displayed as a series of background scenes in the animation environment (Fig. 2c & d i-iii). Next, the 3-D models were imported into the animation software and manually manipulated by the user to “register” the models to their outlines visible in the background images (Fig. 2c & d). Lastly, inter-cube translational and rotational displacements were calculated using this technique. All measurements required for fabricating the grid and solid-body cubes and for measuring the experimental displacements during scanning were performed by a digital caliper (sensitivity = 0.02 mm) (Global Industrial, Port Washington, NY, USA). The details of each step are described below.
Calibrating the MR Imaging Volume: The volume of the MRI coil was calibrated using a custom-built calibration grid (Fig. 2a). Four square Perspex fiber plates (area = 80 mm2; thickness = 2 mm) (Modular-Movement Tray-Set, Games Workshop/NG, UK) were serially stacked with a distance of 30 mm between each plate with three wooden dowels drilled across the plates and glued at all their contact points for stability. Before fixing the dowels, sixteen holes, each 2 mm in diameter, were drilled into each plate in a 4X4 array. Adjacent holes were drilled 20 mm apart from each other. Each hole was fitted with a 2-mm-diameter water bead using a small amount of glue. Three additional beads were embedded into two adjacent plates to define x, y, and z coordinates of the grid (Fig. 2a) . The x- and z-axes were located in the same plane representing the plane of the grid plates, whereas the y-axis extended perpendicular to the plane of grid plates (Fig. 2a). These coordinates were assigned as per the joint coordinate system (JCS) defined by the Standardization and Terminology Committee of the International Society of Biomechanics for studying inter-vertebral motion . To facilitate visualization of the beads in the MR images, the grid was submerged in a 1 % saline solution for 30 s and then air dried for 2 min prior to scanning. The y-axis of the grid was placed along the longer axis of the MRI coil bore (DPA Wrist Coil, Esaote, Genoa, Italy). Four non-contiguous axial 3-mm-thick slices were acquired parallel to and across the grid in a way that each slice image included a single plate with all the 16 beads of a plate in view using a Fast Spin Echo T2 sequence (repetition time [TR] = 7810 ms, time to echo [TE] = 120 ms, field of view [FOV] = 200 x 200, Matrix = 256 × 256; resolution = 0.78 mm; voxel dimension = 1.82 mm3). The four DICOM files were then transferred to the AVIZO software, where all the beads were segmented and images of all segmented individual plates were stored in the TIF format using Photoshop software (Adobe Systems Inc., San Jose, California) for later use in Step 4. Additionally, the surface rendition of the segmented beads representing a composite view of the entire grid was saved as an .OBJ file for digitization in Step 4.
Constructing and Segmenting the Solid-Body Cubes: Two solid-body cubes, with sides measuring ~40 mm, were cut from a wood block (Fig. 2b). Hourglass shaped holes (7-mm base diameters) were drilled through the center of both cubes with a stepped-cone drill. These holes were drilled to create a distinct morphological feature within the cubes and to facilitate the rotoscoping and model-to-scene matching process in a later step. The cubes measured close to the average transverse dimensions of the first lumbar vertebral body in humans, and the hourglass feature simulated the appearance of the vertebral canal in a motion segment [19, 66–68]. Adjacent edges of the cubes were marked with a 20-mm scale with 1.0 mm graduations (Fig. 2b). A neutral position was defined as zero displacement between the mid-lines of the scales. The relative positions between these mid-lines on the scales were manipulated by the user to perform the translation trials with a range of 20 mm in either direction of the neutral position. The opposite sides of the cubes were marked with a protractor to measure inter-cube rotations on both sides of a 0° neutral position at increments of 5° through 90° of rotational displacement. Additionally, 3-D cube models were manually segmented in AVIZO using contiguous high resolution (pixel = 0.78 mm) axial T1 weighted scans (TR = 810 ms, TE = 30 ms, FOV = 200 × 200, Matrix = 256 × 256) (Fig. 2b).
Displacement Trials: The solid-body cubes were immersed in ~1 % saline for 30 s, wiped dry, and positioned within the MR coil. The long axes of the hourglass-shaped holes in both the cubes were placed along the y-axis of the imaging volume and scanned in the neutral position. The single-plane translations and rotations were performed in the z-plane of the imaging volume. The axis for the rotation trials was formed by the x-plane. The cubes were placed and fixed by double-sided tape on a flat foam platform in the coil to avoid shifting during scans. After scanning the neutral position, the platform was pulled out of the coil; and the cubes were re-positioned for the next trail, with the displacements verified by the Vernier caliper before the platform was re-positioned inside the MRI coil (Fig 2b). A gap of ~10 mm was maintained between adjacent edges of the cubes during translation, a distance representing the average dimension of a human lumbar disc space . For the rotation trials, the center of rotation (COR) of the rotating cube was kept 50 mm away from the center of the stationary cube. A high-resolution (0.78 mm) T1 weighted sequence (TR = 810 ms, TE = 30 ms, FOV = 220 × 220, Matrix = 256 × 256, slices = 3, gap = 0, thickness = 5 mm, scan time = ~2 mins/scan) and a fast contrast-enhanced streaming sequence with resolution 0.98 mm (2D hybrid contrast enhanced streaming sequence [2D HYCE S]; thickness = 8 mm, slice = 1, scan time = ~10 s/scan) were used to acquire single-slice images of displaced positions in the mid-sagittal (zy-) plane with the central core of both cubes in view (Fig. 2c). The trials included translations between 0.0–20 mm in 5 mm increments (n = 35 trials; 7 trials/displacement) and rotations ranging between 0o to 20o in 5o increments (n = 30 trials; 6 trials/displacement) on both sides of the neutral position (Table 1). Biplanar translations were scanned in static positions after the cubes were displaced both in the z- and x-planes through a range of 5 mm increments in a 0.0 to 20 mm range (n = 20; 4 trials/displacement). All trials were number-coded and randomly performed in three separate scanning sessions each designated for translational, rotational, and biplanar trials, respectively. For the biplanar trials, additional orthogonal slices were acquired with the central parts of both cubes in view.
Animation of the Imaging Volume and Quantification in MAYA [6, 21, 70]: The MAYA software was used to create the animation environment. The environment essentially represented the calibrated MR imaging volume. The software also provided a “camera-view” for the user to view the cube models and the background scene in the calibrated animation environment (Fig. 2d). The user manually manipulated the 3-D models to match and register them to their outlines seen in the background scenes. The steps for creating the animation environment are as follows:
Creating a MAYA framespec File :
The composite grid .OBJ file created in Step 1 was transferred to MAYA, all the beads were serially numbered according to their actual positions in the grid system, and the centroid points for each segmented bead was calculated by the program. Next, the values of the coordinate points for each bead centroid were calculated in the context of all other beads, representing the entire grid volume. The x, y, and z values of all bead coordinate points were merged together to generate the MAYA ‘framespec’ file to be used for the next step of grid digitization.
Digitizing the Beads:
A MEL-script (MAYA Embedded Language-script) command was run in Matlab. An image of a grid plate previously segmented in AVIZO and stored in a TIF format in Step 1 was opened using the Matlab program, and all beads in the plate were digitized serially by clicking over their central points. Next, the framespec file created in the previous step was loaded into the program to yield the Direct Linear Transformation (DLT) coefficient values for the concerned plate [6, 21]. All four plates were digitized sequentially to generate the DLT coefficient value for each plate image. The program allowed automated corrections for minimization of errors and to contain coefficient values ≤ 1 . This step was repeated for all four plates, and each step yielded a plate coefficient value and a “xyz point” .csv file specific to the concerned digitized plate. The data points of the xyz-point files from all the four plates were collated to generate a common “four-plate xyz point” .csv file for the grid . Next, the four-plate xyz-point and the framespec files were loaded into the Matlab program using the MEL-script. One of the segmented plate-images were opened in the Matlab and re-digitized. The four-plate xyz-point file was loaded into the program, and the MEL-script was re-run to generate a “MayaCam”.csv file that was used to re-create the MR imaging within the animation software and to create the camera-view for the user.
Rotoscoping, Animation, and Quantification:
The animation scene was created using the MayaCam file. After the animation environment was created, the background scene was introduced by importing the TIF format images of the trials into MAYA (Fig. 2c & d). These images were clustered into a series of frames, with each series representing a specific trial type. Next, the two 3-D cube models were imported into MAYA using a scaling factor of 0.1 (from segmentation environment in mm to the animation environment in cm). The models were manipulated with the computer mouse and keyboard functions to achieve maximum geometrically alignment and match between the 3-D model and corresponding image outlines in the background scene. The sharp external boundaries and outlines of the hourglass silhouettes within the solid bodies were utilized to facilitate the model-to-scene match (Fig. 2c & d). This process was called rotoscoping. Once rotoscoping in the first frame of the scene (neutral position in the series) was achieved, an Anatomical Joint Axis (AJX) was assigned to the solid bodies. The image of the background scene was then advanced to the next frame and the rotoscoping repeated; this process was repeated for all remaining trial images. For the biplanar translations, two orthogonal camera views were created to provide background scenes of displacements from two different, the x- and z-plane, perspectives (Fig. 2d). Although the animation software generated solid-body motion data for all rotoscoped image frames in all six degrees of freedom, only applicable single-plane measures were extracted for analysis and reporting. Two sessions (S1 and S2) of rotoscoping and displacement quantification were performed separately for translational, rotational, and biplanar motion by a single observer (NKM). All trials were number coded, and the rater was blinded to the displacement type and the pulse sequence used for the scan. The AJX created in Session 1 was used for rotoscoping in the corresponding Session 2. The approximate time for rotoscoping (Step 4(c)) a series of image frames representing a specific trial type, e.g., a seven-translation series, in this study was ~40 min including matching of the neutral position at the start and extracting displacement data from the series at the end.
Test–retest reliability for the outcomes involving the T1 and 2D HYCE S sequences from the two sessions were determined by coefficient of variation (CV), t-test, intraclass correlation coefficients (ICC) (two-way random effects model with a single measure of reliability), and 95 % limits of agreement (LOA) analyses. Variability between the outcomes from a single displacement quantified in two different sessions was analyzed using CV. For example, if a particular displacement was quantified as 11.5 mm in Session 1 and 10.8 mm in Session 2, the CV was calculated as: Standard Deviation of the two sessions divided by the mean of the two sessions times 100. Thus, for this example, the CV = (0.5/11.15) × 100 = 4.44 %. The sessions were performed at an interval of one week. Additionally, we used dependent sample t-tests to compare the values between testing sessions.
The ICC was calculated using a (2, 1) two-way random effects model with a single measure of reliability computed over the variance observed in the two sessions. A (2, 1) model was chosen as it allows the determination of any existing systematic bias. The statistical software SPSS (SPSS Inc., Version 21.0, Chicago, IL, USA) was used to calculate the ICC. The main objective of the statistical analysis was to ascertain the reliability of this technique. The relative reliability was assessed by calculating the ICC, which assesses the reproducibility of a measurement relative to a sample of repeated measurements. The absolute stability of a measure typically defines the contribution of the main-error component in the observed variance. To fully understand the absolute stability or reliability of a measure, it is essential to understand the contribution of different components of the measurement error [71, 72]. Accordingly, the measurement error was broken into two components. The first component was defined as the systematic bias, and the second was termed the random error. The systematic bias denoted the contribution of any learning effect on the part of the assessor in explaining the between-session variability of the data, whereas the random error explained a biological or mechanical effect [72, 73]. The first step was to generate Bland-Altman plots using the between-session means and differences data. The correlation (R2) between the absolute differences and the means of the between-session values was calculated to determine the spread of the dependent variable. R2 values between 0 and 0.1 represented homoscedasticity, suggesting that there was no correlation between the size of the error and the size of the measured variable. Heteroscedasticity was considered to be present with R2 values > 0.1 and indicated that the degree of error increased with increase in the values of the measured variable along the scale, e.g., the error term increased as the technique attempted to measure larger displacements (translation/rotation) in the experiment [71–74]. Finally, the ratio LOA was calculated for verification of the absolute reliability of the measure using the following equation: ratio LOA = [(SDdiffs/AVGmeans) × 1.96] × 100, where SDdiffs was defined as the Standard Deviation of the difference of scores (Session 1 and Session 2), AVGmeans represented the average of the mean scores (Session 1and Session 2) for each measurement, and the factor 1.96 specified the inclusion of 95 % of observations of the differences in scores. The ratio LOA was interpreted as the highest percentage by which two tests will differ due to measurement error in either the positive or negative direction .
Summary of results
Descriptive statistics and CV and ICC reliability measures obtained from the T1 and 2D HYCE S images for each type of displacement are provided in Table 1. A high degree of between-session reliability was observed for both the T1 and the contrast-enhanced dynamic pulse sequences. Specifically, the average CVs for translation were 4.31 % and 5.26 % for the two pulse sequences, respectively, while the ICCs were 0.99 for both sequences. For rotation measures, the CVs were 3.19 % and 2.44 % for the two pulse sequences with the ICCs being 0.98 and 0.97, respectively. A novel biplanar imaging approach also yielded high reliability, with mean CVs of 3.39 % and 2.66 % noted for translation in the z- and x-plane, respectively, along with ICCs of 0.97 in both planes. Additionally, all but one displacement variables showed homoscedastic relationships with the Bland-Altman’s LOA analysis of the between-session measurements and demonstrated a relatively low degree of systematic bias.
Analysis of the between-session measurements of each of the two sequences, applying a paired sample 2-tailed t-test, did not show any significant differences in means between the T1 (p > 0.98) and the 2D HYCE S (p > 0.84) pulse sequences. The reliability of the measured variables demonstrated high levels of consistency, with the T1 sequence having CVs ranging from 1.1 to 14.63 % and an ICC of 0.99. For the 2D HYCE S sequence, the CVs ranged from 1.6 to 13.5 %, and the ICC was also 0.99. The Bland-Altman plot with 95 % confidence interval (±1.96*standard deviation [SD]) analysis of the between-session data showed that all cases had a test-retest difference within ±1.24 mm (mean/bias = 0.02 mm) and ±1.59 mm (mean/bias = -0.34 mm) for the T1 and the 2D HYCE S sequences, respectively (Fig. 3a). The LOA analysis for translation indicated a relatively low degree of systematic bias in the between-session differences (p = 0.98 and 0.84) and a homoscedastic relationship between the differences and averages of the between-session measurements for both T1 and the 2D HYCE S pulse sequences respectively (R2 = 0.07 and R2 = 0.03) (Fig. 3a). The homoscedasticity indicated that the random errors did not increase with the increase of the measured displacement values. The follow-up ratio LOA analysis demonstrated a systematic bias in the order of 0.02 and -0.34 and random error of ±14.14 and ±13.68 for the T1 and 2D HYCE S sequences, respectively. The ratio LOA analysis for translation suggested that the between-session measurement errors obtained with the technique did not exceed 14.15 % and 13.34 % in either the positive or negative direction with the use of T1 and the 2D HYCE S pulse sequences, respectively.
Analysis of the between-session measurements of the two sequences applying a paired sample 2-tailed t-test did not show any significant mean differences for the T1 (p > 0.94) and the 2D HYCE S (p > 0.96) sequences. The reliability of the measured variables demonstrated high levels of consistency, with the T1 sequence having CVs ranging from 1.2 to 7.6 % and an ICC of 0.98. For the 2D HYCE S sequence, the CVs ranged from 1.05 to 7.6 %, and the ICC was 0.98. The Bland-Altman plot with 95 % confidence interval (±1.96*SD) analysis of the between-session data showed that all cases had a test-retest difference within ±1.27° (mean/bias = -0.14°) and ±0.65° (mean/bias = 0.09°) for the T1 and the 2D HYCE S sequences, respectively. The LOA analysis for the rotation trials indicated a relatively low degree of systematic bias in the between-session differences (p = 0.94 and 0.96) and a homoscedastic relationship between the differences and averages of the between-session measurements for both T1 and the 2D HYCE S pulse sequences respectively (R2 = 0.06 and 0.04) (Fig. 3b). A homoscedastic relationship indicated that the random errors did not increase with the increase of the measured values. The follow-up ratio LOA analysis indicated a systematic bias in the order of -0.15 and -0.09 and random error of ±14.55 and ±20.10 for the T1 and 2D HYCE S sequences, respectively. The ratio LOA analysis for rotation suggested that the between-session measurement errors obtained with the technique did not exceed 14.40 % and 20.01 % in either the positive or negative direction using the T1 and the 2D HYCE S pulse sequences, respectively.
Comparing outcomes between the pulse sequences
Analysis of the difference between the averages of Session 1 and Session 2 translations obtained by T1 and 2D HYCE S pulse sequences did not show any significant results using an independent sample 2-tailed t-test (p = 0.83). The Bland-Altman plot with 95 % confidence interval (±1.96*SD) analysis of the between-sequence data showed a test-retest difference within ±1.50 mm (mean/bias = 0.35 mm) (Fig. 4a). A small heteroscedastic relationship observed in the translation measures indicated that the T1 vs 2D HYCE S between-sequence difference in measured translations increased with assessments of larger magnitudes of translation (R2 = 0.24). The follow-up ratio LOA analysis demonstrated a systematic bias in the order of 0.354 and random error of ±13.41. The ratio LOA analysis for translation suggested that between-sequence measurement errors were within 13.77 % in either the positive or negative direction.
Analysis of the difference between the averages of Session 1 and Session 2 rotations obtained by T1 and 2D HYCE S pulse sequences did not show any significant results using an independent sample 2-tailed t-test (p = 0.98). The Bland-Altman plot with 95 % confidence interval (±1.96*SD) analysis of the between-sequence data showed a test-retest difference within ±0.95° (mean/bias = 0.02°) (Fig. 4a). A homoscedastic relationship observed in the rotation measures indicated that the between-sequence random errors did not increase with assessments of larger magnitudes of translation (R2 = 0.03). The follow-up ratio LOA analysis demonstrated a systematic bias in the order of 0.03 and random error of ±14.28. The LOA ratio analysis for rotation suggested that the between-sequence measurement errors with the T1 and 2D HYCE S pulse sequences were within 14.31 % in either the positive or negative direction.
Analysis of the between-session measurements with a paired sample 2-tailed t-test did not show any significant difference for the z- (p = 0.79) and x (p = 0.73) planes. The reliability of the measured variables demonstrated high levels of consistency, with the z-plane having CVs ranging from 1.51 to 6.15 % and an ICC of 0.97. The CVs for the x-plane measurements ranged from 1.32 to 3.95 %, and the ICC was 0.97. The Bland-Altman plot with the 95 % confidence interval (±1.96*SD) analysis of the biplanar between-session data showed that all displacements had a test-retest difference within ±1.41 mm (mean/bias = -0.04 mm) and ±2.70 mm (mean/bias = 0.12) for the x- and z-planes, respectively (Fig. 4b). The LOA analysis for the biplanar trials demonstrated a relatively low degree of systematic bias in the between-session differences (p = 0.79 and 0.73) and a systematic bias in the order of -0.04 and -0.01 and random error of ±10.87 and ±21.03 for the z- and for the x-planes, respectively. The ratio LOA analysis for translation suggested that the between-session measurement errors obtained with the technique did not exceed 21.03 % and 10.76 % in either the positive or negative direction for the x- and z-planes, respectively. The ratio LOA analysis for the biplanar translation trials showed a homoscedastic relationship between the differences and averages of the between-session measurements for the z-plane (R2 = 0.04), indicating that the random errors did not increase with the increase of the measured values. The x-plane data, however, showed a marginal heteroscedastic relationship (R2 = 0.11), indicating that the random errors did marginally increase with the increase of the measured values (Fig. 4b).
In this study, we describe a novel MRI-based approach that is conceptually similar to some fluoroscopy-based modeling protocols with the major difference of not requiring exposure to ionizing radiation, which has obvious implications for clinical research [6, 18, 19]. While this is the first step in the development of an MRI-based protocol of this nature, our initial work indicates that this technique has promise as we have successfully developed a logical and rational approach to the quantification of motion and have also demonstrated relative and absolute reliability. Below we discuss our findings within the context of the extant literature as well as our future directions.
As stated above, the primary innovation of this work is that it represents an MRI-only, morphology-based modeling technique for tracking solid-body displacements, which is similar to fluoroscopy-based approaches, such as RSA, SR, and XROMM, and static x-ray-based techniques. The scope of application of these fluoroscopy-based techniques is limited due to the ionizing radiation exposure. For instance, obtaining serial measures involving significant radiation exposure over time in research studies requiring oversight by an institutional review board (or other analogous committees charged with approving, monitoring and reviewing biomedical research involving humans) could raise questions about the cost-to-benefit ratio, particularly in light of the Institute of Medicine’s recommendation on avoiding unnecessary medical radiation throughout life [26, 28]. Accordingly, we believe that an MRI-only-based modeling technique for investigating joint kinematics has significant advantages, particularly for the advancement of clinical research.
Available MRI modeling techniques have usually applied multi-slice imaging of the objects of interest to capture the experimental displacements introduced into these objects in the scanning environment. Our study has uniquely attempted a morphology-based single-plane and an orthogonal imaging protocol to quantify experimentally induced displacements in our models. Additionally, we have used a fast-scanning protocol with dynamic contrast-enhanced pulse sequence and compared its outcome to a standard high-resolution T1 imaging. Both these methods have demonstrated high-levels of reliability in quantifying displacements in objects within the MR imaging volume. These findings provide basic proof-of-concept for the notion that a reliable non-ionizing-radiation-based motion quantification technique can potentially be used to characterize a quasi-static visualization of joint kinematics from a single and biplanar approach. The use of dynamic sequence and image processing can be further explored to attempt quantification of joint kinematics in synchronized motion. Additionally, while our single-plane technique does not objectively address detecting out-of-plane motion, inclusion of the orthogonal imaging in the biplanar approach helps manual positioning of the model to match the corresponding out-of-plane shifts of the image silhouettes.
While our initial development results are promising, our study has some limitations. First, we have used static two-dimensional imaging for quantification purposes; and we do not know whether comparable levels of reliability would have been observed if the dynamic pulse sequence were used to scan the solid-bodies in real-time during an un-synchronized motion with subsequent quantification of these images using the technique reported here. The approach we chose was based on technology currently available; to our knowledge, an MRI-compatible device that would permit real-time manipulation of motion is not commercially available, and the custom development of such a device would require significant resources. Second, the current approach required manual segmentation and post-processing, which is very time intensive. Accordingly, we do not know how the use of semi-automatic protocols or automatic iterative segmentation algorithms would have changed our results. Third, while we reported high levels of reliability for a novel biplanar imaging modality (i.e., quantification of motion in two planes, or coupled motion), the orthogonal images obtained for this analysis were not acquired simultaneously (i.e., an image slice was first acquired in one plane and then acquired for the corresponding orthogonal plane) due to the inherent limitation of MR imaging to do so. While this is not necessarily a limitation of the current work, it could pose a limitation for future work that seeks to acquire simultaneous multi-planar images of motion. Lastly, we only assessed reliability and did not assess accuracy. We are currently conducting experiments that will assess accuracy of our technique in a porcine spine model.
In summary, this work provides basic proof-of-concept for a reliable non-ionizing-radiation-based marker-less imaging technique that can potentially be used to quantify quasi-dynamic displacements between joint elements. Additionally, this morphology-based MRI-only technique could be explored further as a tool for real-time joint kinematics analysis.
magnetic resonance imaging
- 2D HYCE S:
Two-Dimensional Hybrid Contrast Enhanced Streaming Sequence
coefficient of variation
digital imaging and communications in medicine
intraclass correlation coefficient
limits of agreement
X-Ray Reconstruction of Moving Morphology
anatomical joint axis
field of view
Cappozzo A et al. Human movement analysis using stereophotogrammetry. Part 1: theoretical background. Gait Posture. 2005;21(2):186–96.
Della Croce U et al. Human movement analysis using stereophotogrammetry. Part 4: assessment of anatomical landmark misplacement and its effects on joint kinematics. Gait Posture. 2005;21(2):226–37.
Reinschmidt C et al. Effect of skin movement on the analysis of skeletal knee joint motion during running. J Biomech. 1997;30(7):729–32.
Cereatti A, Della Croce U, Cappozzo A. Reconstruction of skeletal movement using skin markers: comparative assessment of bone pose estimators. J Neuroeng Rehabil. 2006;3:7.
Chiari L et al. Human movement analysis using stereophotogrammetry. Part 2: instrumental errors. Gait Posture. 2005;21(2):197–211.
Brainerd EL et al. X-ray reconstruction of moving morphology (XROMM): precision, accuracy and applications in comparative biomechanics research. J Exp Zool A Ecol Genet Physiol. 2010;313(5):262–79.
Yamazaki T et al. Improvement of depth position in 2-D/3-D registration of knee implants using single-plane fluoroscopy. IEEE Trans Med Imaging. 2004;23(5):602–12.
Axelsson P, Karlsson BS. Intervertebral mobility in the progressive degenerative process. A radiostereometric analysis. Eur Spine J. 2004;13(6):567–72.
Selvik G. Roentgen stereophotogrammetric analysis. Acta Radiol. 1990;31(2):113–26.
Bojan AJ et al. Three-dimensional bone-implant movements in trochanteric hip fractures: Precision and accuracy of radiostereometric analysis in a phantom model. J Orthop Res. 2015;33(5):705–11.
Bey MJ et al. Accuracy of biplane x-ray imaging combined with model-based tracking for measuring in-vivo patellofemoral joint motion. J Orthop Surg Res. 2008;3:38.
Martin DE et al. Model-based tracking of the hip: implications for novel analyses of hip pathology. J Arthroplasty. 2011;26(1):88–97.
Lai JY et al. A new registration method for three-dimensional knee nearthrosis model using two X-ray images. Comput Methods Biomech Biomed Engin. 2010;13(2):265–78.
Hurschler C et al. Comparison of the model-based and marker-based roentgen stereophotogrammetry methods in a typical clinical setting. J Arthroplasty. 2009;24(4):594–606.
Bey MJ et al. Validation of a new model-based tracking technique for measuring three-dimensional, in vivo glenohumeral joint kinematics. J Biomech Eng. 2006;128(4):604–9.
Dial KP, Goslow GE, Jenkins FA. The Functional-Anatomy of the Shoulder in the European Starling (Sturnus-Vulgaris). J Morphol. 1991;207(3):327–44.
Jenkins Jr FA, Dial KP, Goslow Jr GE. A cineradiographic analysis of bird flight: the wishbone in starlings is a spring. Science. 1988;241(4872):1495–8.
Zhu Z et al. The accuracy and repeatability of an automatic 2D–3D fluoroscopic image-model registration technique for determining shoulder joint kinematics. Med Eng Phys. 2012;34(9):1303–9.
Baier DB, Gatesy SM. Three-dimensional skeletal kinematics of the shoulder girdle and forelimb in walking Alligator. J Anat. 2013;223(5):462–73.
Gatesy SM, David BB. The Origin of the Avian Flight Stroke: A Kinematic and Kinetic Perspective. Paleobiology. 2005;31(3):382–399.
Gatesy SM et al. Scientific rotoscoping: a morphology-based method of 3-D motion analysis and visualization. J Exp Zool A Ecol Genet Physiol. 2010;313(5):244–61.
Anderst WJ, Vaidya R, Tashman S. A technique to measure three-dimensional in vivo rotation of fused and adjacent lumbar vertebrae. Spine J. 2008;8(6):991–7.
Tashman S et al. Dynamic function of the ACL-reconstructed knee during running. Clin Orthop Relat Res. 2007;454:66–73.
Schueler BA. The AAPM/RSNA physics tutorial for residents: general overview of fluoroscopic imaging. Radiographics. 2000;20(4):1115–26.
Latini G et al. Reproductive effects of low-to-moderate medical radiation exposure. Curr Med Chem. 2012;19(36):6171–7.
Kesavachandran CN, Haamann F, Nienhaus A. Radiation exposure of eyes, thyroid gland and hands in orthopaedic staff: a systematic review. Eur J Med Res. 2012;17:28.
Hricak H et al. Managing radiation use in medical imaging: a multifaceted challenge. Radiology. 2011;258(3):889–905.
Mellor FE, Thomas P, Breen A. Moving back: The radiation dose received from lumbar spine quantitative fluoroscopy compared to lumbar spine radiographs with suggestions for dose reduction. Radiography. 2014;20(3):251–7.
Boice Jr JD, Mandel JS, Doody MM. Breast cancer among radiologic technologists. JAMA. 1995;274(5):394–401.
Licciardone JC, Gatchel R, Dagenais S. Assessment and management of back pain. JAMA Intern Med. 2014;174(3):478–9.
Freddolini M, Strike S, Lee RY. Stiffness properties of the trunk in people with low back pain. Hum Mov Sci. 2014;36:70–9.
Hasegewa K et al. Biomechanical evaluation of segmental instability in degenerative lumbar spondylolisthesis. Eur Spine J. 2009;18(4):465–70.
Pearcy MJ, Bogduk N. Instantaneous axes of rotation of the lumbar intervertebral joints. Spine (Phila Pa 1976). 1988;13(9):1033–41.
Gertzbein SD et al. Centrode characteristics of the lumbar spine as a function of segmental instability. Clin Orthop Relat Res. 1986;208:48–51.
Panjabi MM. Clinical spinal instability and low back pain. J Electromyogr Kinesiol. 2003;13(4):371–9.
Dvorak J et al. Clinical validation of functional flexion-extension roentgenograms of the lumbar spine. Spine (Phila Pa 1976). 1991;16(8):943–50.
McGill SM, Cholewicki J. Biomechanical basis for stability: an explanation to enhance clinical utility. J Orthop Sports Phys Ther. 2001;31(2):96–100.
Passias PG et al. Segmental lumbar rotation in patients with discogenic low back pain during functional weight-bearing activities. J Bone Joint Surg Am. 2011;93(1):29–37.
Weiler PJ, King GJ, Gertzbein SD. Analysis of sagittal plane instability of the lumbar spine in vivo. Spine (Phila Pa 1976). 1990;15(12):1300–6.
Yeager MS, Cook DJ, Cheng BC. Reliability of computer-assisted lumbar intervertebral measurements using a novel vertebral motion analysis system. Spine J. 2014;14(2):274–81.
Shin JH et al. Investigation of coupled bending of the lumbar spine during dynamic axial rotation of the body. Eur Spine J. 2013;22(12):2671–7.
Tan Y et al. Kinetic magnetic resonance imaging analysis of lumbar segmental mobility in patients without significant spondylosis. Eur Spine J. 2012;21(12):2673–9.
Yao Q et al. Motion characteristics of the lumbar spinous processes with degenerative disc disease and degenerative spondylolisthesis. Eur Spine J. 2013;22(12):2702–9.
Pearcy MJ. Stereo radiography of lumbar spine motion. Acta Orthop Scand Suppl. 1985;212:1–45.
Pearcy MJ, Hindle RJ. New method for the non-invasive three-dimensional measurement of human back movement. Clin Biomech (Bristol, Avon). 1989;4(2):73–9.
Schneider G, Pearcy MJ, Bogduk N. Abnormal motion in spondylolytic spondylolisthesis. Spine (Phila Pa 1976). 2005;30(10):1159–64.
Yang JS et al. Dynamic Radiographic Results of Different Semi-Rigid Fusion Devices for Degenerative Lumbar Spondylolisthesis: "Dynamic Rod" vs. "Dynamic Screw Head". Turk Neurosurg. 2016;26(2):268–73.
Keorochana G et al. Effect of sagittal alignment on kinematic changes and degree of disc degeneration in the lumbar spine: an analysis using positional MRI. Spine (Phila Pa 1976). 2011;36(11):893–8.
Lee SH, Daffner SD, Wang JC, Davis BC, Alanay A, Kim JS. The change of whole lumbar segmental motion according to the mobility of degenerated disc in the lower lumbar spine: a kinetic MRI study. Eur Spine J. 2014. doi:10.1007/s00586-014-3277-z.
Chung SS et al. Effect of low back posture on the morphology of the spinal canal. Skeletal Radiol. 2000;29(4):217–23.
Rebmann AJ, Sheehan FT. Precise 3D skeletal kinematics using fast phase contrast magnetic resonance imaging. J Magn Reson Imaging. 2003;17(2):206–13.
Behnam AJ, Herzka DA, Sheehan FT. Assessing the accuracy and precision of musculoskeletal motion tracking using cine-PC MRI on a 3.0 T platform. J Biomech. 2011;44(1):193–7.
Ward SR et al. Assessment of patellofemoral relationships using kinematic MRI: Comparison between qualitative and quantitative methods. J Magn Reson Imaging. 2002;16(1):69–74.
Quick HH et al. Real-time MRI of joint movement with TrueFISP. J Magn Reson Imaging. 2002;15(6):710–5.
Shellock FG et al. Kinematic magnetic resonance imaging of the joints: techniques and clinical applications. Magn Reson Q. 1991;7(2):104–35.
Clarke EC et al. A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo: demonstration in the human ankle joint and Achilles tendon. Med Eng Phys. 2015;37(1):93–9.
Fellows RA et al. Repeatability of a novel technique for in vivo measurement of three-dimensional patellar tracking using magnetic resonance imaging. J Magn Reson Imaging. 2005;22(1):145–53.
Fellows RA et al. Magnetic resonance imaging for in vivo assessment of three-dimensional patellar tracking. J Biomech. 2005;38(8):1643–52.
d'Entremont AG et al. Do dynamic-based MR knee kinematics methods produce the same results as static methods? Magn Reson Med. 2013;69(6):1634–44.
Whitehurst DG, et al. Implementing Stratified Primary care Management for low Back Pain: Cost Utility Analysis alongside a Prospective, Population-based, Sequential Comparison Study. Spine (Phila Pa 1976). 2015.
Manchikanti L et al. Epidemiology of low back pain in adults. Neuromodulation. 2014;17 Suppl 2:3–10.
Schaefer C et al. Pain severity and the economic burden of neuropathic pain in the United States: BEAT Neuropathic Pain Observational Study. Clinicoecon Outcomes Res. 2014;6:483–96.
North RB et al. A review of economic factors related to the delivery of health care for chronic low back pain. Neuromodulation. 2014;17 Suppl 2:69–76.
March L et al. Burden of disability due to musculoskeletal (MSK) disorders. Best Pract Res Clin Rheumatol. 2014;28(3):353–66.
Wu G et al. ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion--part I: ankle, hip, and spine. International Society of Biomechanics. J Biomech. 2002;35(4):543–8.
Cheung JP, Shigematsu H, Cheung KM. Verification of measurements of lumbar spinal dimensions in T1- and T2-weighted magnetic resonance imaging sequences. Spine J. 2014;14(8):1476–83.
Ostrofsky KR, Churchill SE. Sex determination by discriminant function analysis of lumbar vertebrae. J Forensic Sci. 2015;60(1):21–8.
Rhoad RC et al. A new in vivo technique for three-dimensional shoulder kinematics analysis. Skeletal Radiol. 1998;27(2):92–7.
Hong CH et al. Measurement of the normal lumbar intervertebral disc space using magnetic resonance imaging. Asian Spine J. 2010;4(1):1–6.
Nyakatura JA, Fischer MS. Three-dimensional kinematic analysis of the pectoral girdle during upside-down locomotion of two-toed sloths (Choloepus didactylus, Linne 1758). Front Zool. 2010;7:21.
Clark BC, Cook SB, Ploutz-Snyder LL. Reliability of techniques to assess human neuromuscular function in vivo. J Electromyogr Kinesiol. 2007;17(1):90–101.
Atkinson G, Nevill AM. Statistical methods for assessing measurement error (reliability) in variables relevant to sports medicine. Sports Med. 1998;26(4):217–38.
Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet. 1986;1(8476):307–10.
Kaya RD, Hoffman RL, Clark BC. Reliability of a modified motor unit number index (MUNIX) technique. J Electromyogr Kinesiol. 2014;24(1):18–24.
All MEL scripts and animation related software tools were used from XROMM resources developed at and made available by the Brown University, RI, USA.
Source of funding
This work was supported by the American Osteopathic Association (AOA) Research Grant # 13-08-682.
The authors declare that they have no competing interests.
NM conceived the study idea and design, performed image acquisition, data analysis, drafting and revision of manuscript. SM supported in designing the study, development of technical strategy, draft of manuscript. JC provided inputs on data analysis and drafting of manuscript. SW helped with the study design and to draft and revise the manuscript. JT assisted with study design and in editing the manuscript. BC contributed with data analysis, major edits on original manuscript and the revision. All authors read and approved the revisions and the final version of the manuscript.
About this article
Cite this article
Mahato, N.K., Montuelle, S., Cotton, J. et al. Development of a morphology-based modeling technique for tracking solid-body displacements: examining the reliability of a potential MRI-only approach for joint kinematics assessment. BMC Med Imaging 16, 38 (2016). https://doi.org/10.1186/s12880-016-0140-1
- Scientific rotoscoping
- Dynamic sequence
- Back pain