Multi-Modal Data Correspondence for the 4D Analysis of the Spine with Adolescent Idiopathic Scoliosis

Adolescent idiopathic scoliosis is a three-dimensional spinal deformity that evolves during adolescence. Combined with static 3D X-ray acquisitions, novel approaches using motion capture allow for the analysis of the patient dynamics. However, as of today, they cannot provide an internal analysis of the spine in motion. In this study, we investigated the use of personalized kinematic avatars, created with observations of the outer (skin) and internal shape (3D spine) to infer the actual anatomic dynamics of the spine when driven by motion capture markers. Towards that end, we propose an approach to create a subject-specific digital twin from multi-modal data, namely, a surface scan of the back of the patient and a reconstruction of the 3D spine (EOS). We use radio-opaque markers to register the inner and outer observations. With respect to the previous work, our method does not rely on a precise palpation for the placement of the markers. We present the preliminary results on two cases, for which we acquired a second biplanar X-ray in a bending position. Our model can infer the spine motion from mocap markers with an accuracy below 1 cm on each anatomical axis and near 5 degrees in orientations.


Introduction
Adolescent idiopathic scoliosis (AIS) is a 3D spinal disorder affecting 1-3% of the population and that evolves during the period of growth [1]. There is currently a lack of knowledge about its etiopathogenesis. Several factors are investigated, such as genetics, hormones or biomechanical disorders [2]. Scoliosis is usually diagnosed and monitored from 2D X-rays [3], but the recent literature shows that the new static 3D characterization methods, such as the severity-index [4] computed from different 3D descriptors of the deformities (torsion index, Cobb angle, axial rotations, ...), provide valuable insights into the progression of scoliosis. Despite the advances in the static characterization, the dynamic behavior of the condition remains poorly studied in the scientific literature. However, motion capture technology provides a promising approach to analyze the three-dimensional movement of patients with scoliosis.
Motion capture provides a superficial analysis of the patient's motion but does not directly capture the actual dynamics of the spine inside the body. Several works [5][6][7] have investigated how to describe the spinal alignments in motion by acquiring 3D trajectories with mocap markers positioned on the palpable spinous process of the vertebrae. These methods, however, are highly sensitive to the palpation task and underestimate the coronal curvatures measurements [8] without numerical correction regarding the true anatomy [6]. In addition, these methods cannot provide a detailed analysis in rotations at each vertebra level.
A complementary approach is the use of a kinematic model in charge of the reconstruction of the vertebra dynamics according to external constraints. Rigid-body modeling is the usual approach, in which the bones are defined as a set of rigid-bodies articulated with joints that allow different degrees-of-freedom (DOFs) in rotation and/or translation. Rigid-body models have already been used in the scope of static treatment correction prediction and analysis [9][10][11]. Overbergh et al. [7] proposed and validated a workflow allowing the reconstruction of a subject-specific model that can be used towards markerbased motion capture analysis. They leveraged the recent advances in low-dose biplanar radiography provided by the EOS Imaging System (EOS Imaging, Paris, France) to validate their kinematic predictions. However, their method was evaluated for adult spinal deformities and requires a sensitive and time-consuming palpation task which cannot be easily implemented on young AIS patients.
The aim of this study is to build a subject-specific kinematic model of patients with AIS that captures both their internal and external specificities. The digital twin can then be driven with mocap markers on the back of the patient, which yields the actual kinematics of the spine. We present a preliminary validation of the predictions against 3D spine reconstructions of two patients who performed a lateral bending.

Collected Data
Our dataset consists of 8 patients aged between 8 and 16 years with AIS (Cobb angle range: 14-68 • ) and without any treatment history. They were included in our study following the IRB CPP Ile de France 2 on 20 July 2020 : n • ID RCB: 2020-A01071-38. All parents and patients received an information letter. Two data modalities were collected in our institution for each subject: • A biplanar X-ray of their trunk made with an EOS imaging system • A surface scan of the back using an Occipital Structure Sensor Mark II (XRPro, LLC, Saratov) During all acquisitions, the patients wore a set of 18 radio-opaque markers positioned in the neighborhood of 5 vertebrae: T1, T4, T7, T10 and L03. The marker placement does not require a precise palpation. The markers can be identified and located in the X-ray images and the surface skin mesh (Figures 1 and 2).  To validate the kinematic predictions of the model, X-rays of two voluntary patients, with Cobb angles of 14 • and 29 • , were captured in left and right lateral bending. These images will be used to quantify the inference of the spine motion inside the body in different poses.

Data Processing
Some manual steps are needed to annotate the images. The first step is to make the semi-automated reconstruction of the 3D geometries of the spine from the biplanar radiographs with the SterEOS software 1.8.99.21R (copyright 2015 EOS Imaging) [12] (EOS Imaging, Paris, France). Then, the radio-opaque markers, visible in the images, are located in 2D and their 3D position is computed using the calibration information available in the DICOMs metadata [13].
From the surface skin mesh, the marker locations are identified. As the surface scan and the EOS reconstruction are defined in a different global frame, we use a rigid registration to bring them in the same frame. Namely, we compute the translation and rotation that minimizes the distance between the 3D markers identified in both modalities. Let us note that as the pose might be slightly different in both acquisitions, we use RANSAC [14] to filter out markers whose positions might have significantly changed with the patients pose. This registration provides a first association between the surface scan and the spine 3D reconstruction.

The Subject-Specific Kinematic Model
To create a 4D numerical twin of the patient, we leverage the Anatoscope technology based on "Anatomy transfer" [10,15,16]. This method deforms an initial anatomic model to capture the internal and external shapes of the patient with rigid and elastic registration processes. The resulting avatar can then be used for biomedical simulations, namely, the parameters of the model can be optimized so that the skin of the model matches the mocap markers while enforcing biomechanical constraints on the spine behavior.
Our kinematic model has N x = 18 articulated rigid-bodies x i defined in positions p i ∈ R 3 and rotations R i ∈ R 3×3 , corresponding to each i thoracolumbar vertebra and the pelvis. In addition each bone is associated with a set of shape parameters s that modify the model geometries. The bones are connected by K = 17 joints of 6 degrees-of-freedoms (DOFs) in translations and rotations, as defined by Ignasiak et al., 2016 [17]. We followed the modification applied by Koutras et al., 2021 [18], allowing for a symmetrical definition of the vertebral motion ( Table 1). The position of the joints is defined in the middle of the segment drawn between the two adjacent endplate centroids. Their orientation is defined accordingly to the inferior rigid body. The model also has a skin surface that is rigged by the articulated rigid bodies. The first step of the registration process is to change the pose x and the shape s parameters of the model so that the models' spine fits the 3D reconstructed spine and the models' skin surface captures the surface scan. As the shape parameters s only capture several deformation models, the obtained geometries do not precisely match the patientspecific geometries. Thus, in a second step, we refine the geometries of the model's vertebrae and skin to match the patient's observations (EOS geometries and surface scan). Once the patient surfaces are captured by the model, the shape parameters are fixed. Let us note that the patients pose during the surface scan is slightly different than the pose in the radiography, thus, the location of the markers m A on the resulting model differs from the marker locations m D on the X-rays. To fix this issue, we transfer the N m markers positions m S ∈ R 3×N m located on the surface scan onto the model skin mesh. We identify the closest mesh face of the model to a marker and define the marker location on the model mesh using the barycentric coordinates of the face vertices. Then, we use a temporary set of pose parameters x that are optimized so that the model markers m A (x ) match the ones in the X-rays m D . This effectively changes the model skin surface to match the pose of the patient in the EOS device.
The marker-based optimization is computed aŝ where the energy E model is a regularization term enforcing anatomic constraints on the joints of the model. The resulting model skin surface matches the pose of the back surface during the X-ray acquisitions. Thus, we create a synced skin and spine model by disregarding the temporal parameters x and associating the current optimized skin to the original model parameters x obtained during the first registration process. As a result, the anatomical model is a numerical twin of the patient, including the skeleton and the skin rigged with common model parameters. The association of the skin and spine is performed on the pose observed during the X-ray acquisitions with the help of the radio-opaque markers. Given a new set of markers, the model parameters can be optimized using Equation (1) and obtain the skin model matching the input markers, as well as a prediction of the spine geometry inside the body.

Accuracy of the Anatomical Model
To evaluate the quality of our model, we compute several metrics related to the accuracy in shape, positions and orientations of the vertebrae. The ground-truth measurements of the vertebrae are based on the 3D annotations provided by the SterEOS software. We compare our model meshes with the reconstructed ones.
The vertebrae location is defined with the center of mass of the vertebra mesh. The euclidean distance between the corresponding ground truth and model meshes is then computed with their mean absolute difference on each anatomical axis. The differences in orientations are given by the 2D projection of a vertebra orientation vector, computed according to the recommendations of the International Society of Biomechanics (ISB) [19] on a given anatomical plane. The resulting angle between the SterEOS measurement and our model is then measured for each vertebra on each plane. To assess the quality of the model vertebrae geometry, we compute the absolute mean and standard deviation of the point-tosurface distances between the model geometries and the EOS 3D reconstructed spines.
As the body surface of the back during the radiograph acquisitions is not available, we evaluate the model fit to radiographs by computing the 3D euclidean distances between the radio-opaque markers on the model and in the X-ray. We also quantify the contribution of the model skin correction step used to reconstruct the pose of the back during the X-ray acquisition.

Validation of the Kinematic Predictions
We validate the motion of the spine inside the body predicted by our model as follows. X-rays from two voluntary subjects with AIS were acquired in different poses in the EOS imaging system. The standard pose, standing with hands on the cheeks, was used to create the digital twin ( Figure 2a). Then, two other poses, right and left lateral bending, were also acquired ( Figure 2b). From the resulting images, the marker positions were identified and located and the 3D spine model reconstructed. The markers of these poses are used to drive the model, and the obtained 3D spine is compared to the reconstructed one. Specifically, given a set of 3D markers in the X-rays m D , the model parameters x can be optimized so that the model markers m A best match the input markers m D by optimizing Equation (1).
Let us note that the 3D reconstruction of the vertebrae in lateral bending is not straightforward due to the overlapping of the different bones on the profile X-ray (Figure 2b right). Thus, the vertebra details needed for the reconstruction are difficult to extract: the computed geometries of the vertebrae do not accurately match the image and do not match the geometries obtained in the standing pose. To overcome this issue, we rigidly registered the model vertebra geometries obtained in the standing position to the reconstructions in bending position. This step computes the optimal rigid transformation of each vertebra to minimize the projected distances of the model geometries in standing and those obtained in bending.
With this procedure, we obtained the 3D rigid location and orientation of each vertebra from the bending images. Thus, we compute the orientation errors with respect to the predictions by using the intrinsic Euler angles defined by the ISB XYZ sequence (coronal, axial, sagittal). The accuracy in position is given as described in Section 2.4.

Results
We evaluated the created digital twin in two settings. We first quantified the capability of the model to capture the data in the standing position, and then we evaluated the precision of the model in the bending position. For each position, we assess the external accuracy, i.e., how well does the model fit the 3D markers, as well as the internal accuracy, i.e., how well does the model capture the shape and position of the 3D spine inside the body? We tested our method to correct the pose of the model skin according to the 3D positions of the markers in the X-rays. In Figure 3, it appears that we were able to have a significant gain in accuracy in the surface reconstruction of the model, reflected by the position of the markers, with the correction step. The average distance error decreases from 9.86 mm (std: 8.10 mm) to 4.47 mm (std: 2.70 mm).

Internal Accuracy
The measurements in positions are made by computing the center of mass of each mesh. The orientations of each thoracic and lumbar vertebra (T01-L05) are produced according to the ISB recommendations. The shape accuracy is given for each vertebra by computing the mean of the absolute point-to-surface distances. The results are detailed for each vertebra on Table 2. The error in positioning T01 is due to a model registration error on a unique patient whose vertebra has a particular shape. Despite this result, the model performs well in capturing the morphological specificities of the patient's spine.

Accuracy of the Subject-Specific Model in Bending
From the X-ray images, the radio-opaque marker positions were manually extracted and their 3D location triangulated. Their positions serve as inputs for the inverse-kinematics problem.
For this experiment, we removed the two most lateral markers at the T04 level, as these markers are subject to the movement of the scapula, which is not included in our model. We can notice that their positions have not correctly fit their target according to the anatomical model for these two subjects with an average 3D error of 6.68 mm (std: 1.25 mm).
Let us note that for one subject, the right bending pose resulted in most of the markers being out of the X-ray frontal plane view. Thus, we do not report the metrics on this case.

External Accuracy
After the marker optimization, the model markers reached their corresponding targets with an average distance error of 4.66 mm (std: 2.14 mm).

Internal Accuracy
We evaluated the accuracy of our predictions in positions by comparing the center of mass for each vertebra. The errors in rotations are given by comparing the intrinsic Euler angles between the ground-truth and predictions according to the XYZ order given by the ISB recommendations [19]. The results are presented for each vertebra in Table 3. Table 3. Accuracy of the predictions (MAE) of the vertebra positions and orientations according to the marker positions in lateral bending.

Positions (mm)
Orientations ( The 3D accuracy of the predicted positions is close to the actual cm on average with 1.07 cm (std: 0.42 cm). It can be noted that the predictions are affected by a global lateral shift in the side of the movement (Figure 4) highlighted on the mediolateral axis (MAE 8.49 mm, std: 5.13 mm). Despite this error, the orientation in the corresponding plane (coronal) is closer to the expectations (4.57 • , std: 4.53 • , Table 3). The greater error in the L05 transverse rotation can be due to a lack in superficial constraints (i.e., markers) in this region.

Discussion
In this study, we presented a semi-automatic workflow that allows the creation of a 4D numerical avatar of patients with AIS that reflects their inner and external anatomy. The inputs were captured using safe and low-dose imaging methods and do not need the sensitive and time-consuming task of palpation that cannot be applied on young subjects in daily clinical usage. An anatomical avatar is deformed in order to capture the internal (spine) and external (skin) specificities of the patient obtained from the different modalities. The model was able to capture the vertebra geometries with a mean error below the mm (0.29 mm, std: 0.28 mm), allowing us to compute several descriptors of the vertebra positions and orientations automatically. One major challenge is to recover the external shape of the patient during the X-ray acquisitions, as the patient's pose is necessarily different from the surface scan acquired separately. A solution can be found in the introduction of 3D sensors during the X-ray acquisition [20]. We proposed a method that leverages the radio-opaque markers, visible in both the X-ray images and surface scan, to correct the pose of the model's back. This additional step in our workflow allows us to capture the change in pose of the patient in standing and to increase the correspondence with the markers in the X-rays from 9.86 mm (8.10 mm) on average to 4.47 mm (std: 2.70 mm). However, we notice the difficulty for our model to fit the lateral markers, particularly on the upper part of the body. This can be explained by the fact that we are not modeling the shoulder girdle, and some markers, positioned near T04, for instance, are placed on the scapula.
In a preliminary study, we validate the predictions of the kinematic model with secondary X-rays of two voluntary patients, who were asked to make lateral bendings (left and right, Figure 2). One capture was rejected due to the low visibility of the markers. The marker positions were used as inputs in the simulator, and the predicted positions and orientations of the vertebrae were compared to the 3D X-ray reconstructions.
The creation of a thoracolumbar kinematic model of patients with the scoliosis condition was also investigated by Overbergh et al., 2020 [7]. They were able to make a subject-specific model of their patients comprising bone geometries and a set of superficial markers. Then, they compared the kinematic predictions of the models against secondary X-rays of the subjects in different poses. However, this workflow was designed towards adult spinal deformity analysis and did not integrate the external surface of the back. We proposed a method that leverages the 3D acquisition of the back to avoid the precise and time-consuming task of palpation. Thus, this protocol can be used by less experienced medical staff.
Our predictions are in the range of values of those reported in [7] in positions and orientations, except for the mediolateral axis measurements. In this case, the predicted spine was affected by a global lateral shift on the side of the bending. The orientations of the vertebrae were in the range of acceptable values with an average accuracy below 6 • on the intrinsic Euler rotations. The reconstructions in the lumbar part of the spine showed a more important error, particularly on the axial rotations. This is also the most flexible part of the spine where the rotations can be underestimated by our kinematic model. Some limitations can be noted such as the joint stiffness that is defined using the Ignasiak et al., 2016 [17] and Koutras et al., 2021 [18] values obtained from healthy adult observations. The optimization processes can be planned in order to estimate more specific kinematic parameters such as the joint stiffness or the marker constraint. The addition of radio-opaque markers at the hips would allow us to improve the model predictions in the lower part of the spine and to validate predictions with the pelvis. Furthermore, the recruitment of more patients, with and without AIS, would provide a more comprehensive overview of the model performances.

Conclusions
Combined with 3D X-ray imaging, the characterization of the spine in motion would provide valuable insights about the scoliosis condition of the patient. This analysis is challenging since there is no non-invasive method to capture the vertebral motion in vivo in the young patient. In this study we investigated the use of a subject-specific kinematic model obtained from low-dose biplanar X-rays, a surface scan and a set of radio-opaque markers. As a preliminary result, we evaluated the kinematic behavior of the model against secondary biplanar X-rays of two patients in lateral bending. We show that the kinematic predictions are close to the radiograph observations with an accuracy near 1 cm in 3D position and 5 • in orientation.
Future work can be considered, such as the addition of markers or the optimization of the physical model constraints. Leveraging bigger cohorts of patients in motion will further allow us to better characterize the individual physical properties of the patients. Ultimately, a comprehensive understanding of scoliosis using 4D avatars will lead to improved scoliosis classification, its diagnosis and treatment with biomedical simulation.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the local committee CPP Ile de France 2 on the 07/20/2020, num. ID RCB 2020-A01071-38.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: Not applicable.