Four-Dimensional CT Analysis Using Sequential 3D-3D Registration

Summary

We analyzed joint kinematics from four-dimensional computed tomography data. The sequential 3D-3D registration method semiautomatically provides the kinematics of the moving bone with respect to the subject bone from four-dimensional computed tomography data.

Abstract

Four-dimensional computed tomography (4DCT) provides a series of volume data and visualizes joint motions. However, numerical analysis of 4DCT data remains difficult because segmentation in all volumetric frames is time-consuming. We aimed to analyze joint kinematics using a sequential 3D-3D registration technique to provide the kinematics of the moving bone with respect to the fixed bone semiautomatically using 4DCT DICOM data and existing software. Surface data of the source bones are reconstructed from 3DCT. The trimmed surface data are respectively matched with surface data from the first frame in 4DCT. These trimmed surfaces are sequentially matched until the last frame. These processes provide positional information for target bones in all frames of the 4DCT. Once the coordinate systems of the target bones are decided, translation and rotation angles between any two bones can be calculated. This 4DCT analysis offers advantages in kinematic analyses of complex structures such as carpal or tarsal bones. However, fast or large-scale motions cannot be traced because of motion artifacts.

Introduction

Joint kinematics have been described using a number of methodologies, such as motion capture sensors, 2D-3D registration, and cadaveric studies. Each method has specific advantages and disadvantages. For example, motion capture sensors can measure fast, large-scale motions using infrared cameras with or without sensors on the subject^1,2. However, these methods measure skin motion to infer joint kinematics, and therefore contain skin motion errors³.

Cadaveric studies have been used to evaluate ranges of motion, instability, and contact areas^4,5,6. This approach can measure small changes in small joints using CT or optical sensors attached directly to bone using pins or screws. Cadaveric models can mainly evaluate passive motions, although multiple actuators have been used to apply external forces to tendons to simulate dynamic motion⁷. Active joint motion can be measured by 2D-3D registration techniques, matching 3DCT images to 2D fluoroscopy images. Although the accuracy of the registration process remains controversial, the reported accuracy is generally high enough for large joint kinematics^8,9. However, this method cannot be applied to small bones or multiple bones in narrow spaces.

In contrast, 4DCT is a dynamic CT method that obtains a series of volumetric data. Active joint motions can be analyzed using this approach¹⁰. This technology provides precise 3D positional data of all substances inside the CT gantry. The 3D joint motions are clearly visualized in a viewer. However, describing joint kinematics from such a series of volume data is still difficult, because all the bones are moving and no landmarks can be traced during the active motions in vivo.

We developed a method for 4DCT analysis that provides the in vivo joint kinematics of the whole bones around the joint during active motions. The aim of this article is to present our method, the sequential 3D-3D registration technique for 4DCT analysis, and show representative results obtained using this method.

Protocol

All methods described here have been approved by the Institutional Review Board of Keio University School of Medicine.

NOTE: Joint kinematics are measured by reconstructing the motion of a moving bone around a fixed bone. For knee joint kinematics, the femur is defined as the fixed bone and the tibia is defined as the moving bone.

1. CT imaging protocol

1. Set up the CT machine. Acquire CT examinations with a 320-detector-row CT system to allow for multiple phases of 3D volume data with 160 mm craniocaudal coverage. For example, in the analysis of knee kinematics, the image acquisition consists of 51 volume scans with a rotation time of 0.275 s, and all images are reconstructed using half reconstruction, so that the temporal resolution is approximately 0.16 s.
2. Use the following scanning parameters: peak tube voltage = 100 kVp; tube current = 40 mA; scan coverage = 160 mm; matrix size= 512 x 512 pixels; and reconstruction section thickness and section interval = 0.5 mm.
3. Place the target joint of the subject inside the CT gantry in the starting position of the 4DCT exam (Figure 1).
4. Before the CT exam, rehearse movements of the joint from the start position to the end position within the required examination time. Ask the subject to move the joint during the 10.275 s scan time and obtain a series of volume data. Store the sequential volume data in DICOM format.
5. Perform static 3DCT of all the target bones and store the data in DICOM format.

2. Surface reconstruction

1. Perform semiautomatic segmentation of 3DCT data (Figure 2A).
   1. Load CT DICOM data by selecting all DICOM files of the static 3DCT data.
   2. Open the label field by clicking Edit New Label Field and check which threshold CT attenuation value is appropriate to extract cortical bone from the source bone. Select materials with CT attenuation values above the threshold. For example, the bone cortex threshold for a young subject is set as 250. Check the label for bone cortex selection and manually modify the demarcation using an editing tool for consistency with the shape of the bone.
   3. Generate the surface data (triangle meshes) from the labeled bone cortex position data (point cloud in the software). Store the surface data by exporting data in Standard Triangulated Language (STL) format.
   4. Click Generate Surface|Apply on the label of the cortical bone. Click File|Export Data As|STL Binary Little Endian to save the surface data in STL format.
2. Perform automatic segmentation of 4DCT volume data (Figure 2B).
   NOTE: Each frame of the DICOM data includes the distribution of the CT attenuation values in the CT gantry.
   1. Set the threshold of the bone cortex as in static CT, and extract geometric data showing CT attenuation values above the threshold from all 51 frames of the 4DCT data using the DICOM reading module in the programming software. Adjust the threshold according to the bone density of the source bone. For example, for the osteoporotic bone, set the threshold lower.
   2. Translate all positional data that have already been obtained in the previous step into a format that can be interpreted by image processing software (e.g., Avizo). In the image processing software, reconstruct all surface data of the point cloud with higher CT attenuation values than the threshold for all 4DCT frames using a batch processing script. The image processing software contains the function to read the script and export the surface data from the DICOM series data automatically. The batch script is shown in the Supplemental Coding File.

3. Image registration

NOTE: In this step, reconstruct the motions of the moving bone with respect to the fixed bone from the raw 4DCT DICOM data.

1. Perform surface registration from static 3DCT to the first frame of the 4DCT.
   1. Trim the bones in a static 3DCT into partial segment data that are included in all frames of 4DCT for use with the iterative closest point (ICP) algorithm¹¹ in the 3D mesh editing software using the Selecting Face function (Figure 3A) by referring 4DCT movie data. The surface data from 4DCT are only partial segments that are included in each volume imagebecause surface registration requires that one surface data point is included in another surface.
   2. Pick three landmarks in the fixed and moving bones that can be easily identified from the trimmed 3DCT surface and the surface data of the first frame of 4DCT in the 3D mesh editing software using the PickPoints function (Figure 3B).
   3. Match the partial fixed and moving bones roughly on the first frame of the 4DCT surface data (Figure 3C) according to the picked landmarks in 3.1.2. Next, perform surface registration using the ICP algorithm¹¹ using the open source software (e.g., VTK).
      NOTE: This process provides homogeneous transformation matrices of the fixed and moving bones from the static 3DCT to the first frame of 4DCT (Figure 3D). These matrices are 4 x 4 matrices consisting of rotation and translation, as shown in Figure 4. The transformation matrix causing the reverse action can also be calculated.
2. Perform sequential surface registration (Figure 5).
   1. Match the partial surfaces of the fixed and moving bone in the first 4DCT frame onto the surface data of the second frame. Next, match the partial surfaces of the i^th frame onto the (i + 1)^th frame of 4DCT sequentially. Repeat this process until the last frame of the 4DCT by programming with use of the ICP module in the open source software.
3. Calculate transformation matrices from the static 3DCT to all frames in 4DCT according to the results of 3.1 and 3.2.
4. Reconstruct moving bone motion with respect to the fixed bone (Figure 6).
   1. Reconstruct the kinematics of the moving bone with respect to the fixed bone from the matrices that represent the transformation from the static 3DCT to each 4DCT frame. Define the coordinate systems of the fixed and moving bones when the rotation parameters are measured (e.g., flexion angle or rotation angle calculated by the Euler/Cardan angle)^12,13,14.

Results

We describe the motion of the tibia during knee extension. The knee joint was positioned in the CT gantry. A triangle pillow was used to support the femur at the starting position. The knee was extended to a straight position over the course of 10 s. Radiation exposure was measured. In addition to 4DCT, static 3DCT of the whole femur, tibia, and patella was performed. Surface data of the whole femur and tibia were reconstructed. The threshold for HU numbers of the bone cortex was set as 2...

Discussion

Our method allows visualization and quantification of the motions of whole bones and provides numerical positional data of the moving bone with respect to the fixed bone from 4DCT data. Many tools have been suggested for measuring joint kinematics. Motion skin markers can analyze total body motions over a long time. However, this method contains skin motion errors³. Joint kinematics should be estimated from the motion of adjacent bones. The 2D-3D registration method uses fluoroscopy and infers 3D ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
4DCT scanner	Canon medical systems (Tochigi, Japan)	N/A	4DCT scan, Static 3DCT scan
AVIZO(9.3.0)*	Thermo Fisher Scientific (OR, USA)		Image processing software. Surface reconstruction from CT DICOM data and point cloud data. * Ryan, T. M. & Walker, A. Trabecular bone structure in the humeral and femoral heads of anthropoid primates. Anat Rec (Hoboken). 293 (4), 719-729, doi:10.1002/ar.21139, (2010).
Meshlab**	ISTI (Pisa, Italy)	N/A	Surface trimming and landmark picking ** MeshLab: an Open-Source Mesh Processing Tool. Sixth Eurographics Italian Chapter Conference, page 129-136, 2008. P. Cignoni, M. Callieri, M. Corsini, M. Dellepiane, F. Ganovelli, G. Ranzuglia
VTK(6.3.0)***	Kitware (New York, USA)	N/A	Iterative Closest Points algorithm. Used in python language programming. *** https://vtk.org
Python(3.6.1)	Python Software Foundation	N/A	DICOM file processing to extract the point cloud from the bone cortex ('dicom.py' module). Calculation of the rotation matrices. (Numpy module) Sequential image regestration using ICP algorithm