Postoperative 3D spine reconstruction by navigating partitioning manifolds

Purpose: The postoperative evaluation of scoliosis patients undergoing corrective treatment is an important task to assess the strategy of the spinal surgery. Using accurate 3D geometric models of the patient’s spine is essential to measure longitudinal changes in the patient’s anatomy. On the other hand, reconstructing the spine in 3D from postoperative radiographs is a challenging problem due to the presence of instrumentation (metallic rods and screws) occluding vertebrae on the spine. Methods: This paper describes the reconstruction problem by searching for the optimal model within a manifold space of articulated spines learned from a training dataset of pathological cases who underwent surgery. The manifold structure is implemented based on a multilevel manifold ensemble to structure the data, incorporating connections between nodes within a single manifold, in addition to connections between di ff erent multilevel manifolds, representing subregions with similar characteristics. Results: The reconstruction pipeline was evaluated on x-ray datasets from both preoperative patients and patients with spinal surgery. By comparing the method to ground-truth models, a 3D reconstruction accuracy of 2 . 24 ± 0 . 90 mm was obtained from 30 postoperative scoliotic patients, while handling patients with highly deformed spines. Conclusions: This paper illustrates how this manifold model can accurately identify similar spine models by navigating in the low-dimensional space, as well as computing nonlinear charts within local neighborhoods of the embedded space during the testing phase. This technique allows postoperative follow-ups of spinal surgery using personalized 3D spine models and assess surgical strategies for spinal deformities. C


INTRODUCTION
Structural spinal deformities like adolescent idiopathic scoliosis (AIS) are complex deformations of the trunk described in 3D, including deviations of the spine in the sagittal plane in addition to asymmetric vertebral deformation. Magnetic resonance imaging (MRI) is very attractive because it is noninvasive for the patient, but it is unfortunately not suitable for a postoperative 3D evaluation due to the ringing artifacts caused by the surgical implants. However, computerized tomography (CT) is a more suitable imaging modality than MR imaging in terms of 3D reconstruction of bony structures, but CT-scans expose children to unacceptable doses of ionizing radiations in order to reconstruct the entire spine geometry (all thoracic and lumbar vertebrae). More importantly, both of these modalities cannot be done with the patient standing. Due to these limitations, biplanar radiography remains the most frequently used modality to assess AIS in 3D as it allows the acquisition of data in the natural posture while exposing the patient to a low dose of radiation. 1 Various approaches were proposed for biplanar x-ray spine reconstruction. A first set of methods focused on super-vised techniques requiring manual interventions. A priori representations of morphologic descriptors obtained from geometrical vertebrae helped to model the overall appearance of the spine, which would be subsequently adjusted from shape projections. 2 The method was improved with the introduction of a set of rule-based inferences to yield an accurate estimate of the vertebra's orientation and 3D location. 3 In previous work, 4 the authors presented a similar approach using parametric models that are obtained from the x-rays vertebral midline. Another set of approaches represented the spine as a set of intervertebral articulations, either to infer the overall geometry of the spinal column 5 or for the estimation of the geometrical shape based on second-order cone programming. 6 In previous work, 7 a hybrid statistical and data feature technique was proposed to model the spine based on high-level primitives. On the other hand, all of these above referenced techniques are dedicated to model the preoperative spine from biplanar xrays, where anatomical structures such as the vertebral bodies are visible without any occlusion from spinal instrumentation.
Surgical treatment involves straightening the anatomical curves using preshaped metallic rods that are fixed with pedicle screws into the spinal body. Once the patient's spine has been instrumented, it is important to enable the longitudinal evaluation of the deformity progression in postoperative examinations. Using 3D models of the spinal column from postoperative radiographs, this helps clinicians to measure changes in the spine's curvature. However little efforts have been placed to reconstruct in 3D spines with instrumentation, even though the analysis of postoperative models is important to evaluate the adopted surgical strategy. A hierarchical statistical model was proposed to achieve this objective. 8 In this later work, a training procedure would partition the observations into withingroup and a between-group structures, which is subsequently incorporated into the 3D reconstruction process. The modeling is described as a second-order linear programming problem that can be solved in an efficient manner. This requires a user to identify initial landmarks on the images and is based on a linear sampling of the underlying pathological variations in order to generate the model.
These past linear statistical models differ from manifold learning paradigms which assumes that high-dimensional data can be embedded into a low-dimensional domain. Unfortunately, traditional methods such as locally linear embedding (LLE) or ISOMAP have difficulty in capturing closed distributions since they seek to "unwrap" manifolds onto a linear plane, which ultimately corrupts the intrinsic topology of the low-dimensional domain. Indeed, methods that globally map higher-dimensional data to an intrinsic space are often limited to open manifolds, 9 which are dedicated to model noncyclical data distributions. 10 Techniques based on higherdimensional mappings with closed manifold preservation or selecting subsections of the manifolds to embed often fail for closed manifolds as shown by Pitelis et al. 11 Conversely, methods anchored on graphical models attempt to maintain the overall manifold appearance which preserves the continuity of the embedding. 11 Learning piecewise linear models of a manifold based on a collection of overlapping linear charts was shown to outperform standard approaches. 12 This technique differs to manifold charting that try to preserve the overall topology of the manifold charts, and therefore does not depend on sectioning the manifold into subsections or mapping the data into a higher-dimensional space. This feature can improve the reliability of the inference process when confronted to a new instance from a manifold of lower dimensionality that was embedded in the higher-dimensional domain.
Using these recent advances in manifold charting, we propose a new biplanar 3D reconstruction technique of in-strumented spines, where the shape reconstruction task is formulated as minimizing an energy term from a manifold embedding of instrumented spine models. The main advantage of the minimizing cost function lies in its tolerance toward spine instrumentation, yielding valid anatomical configurations based on a set of prior representations. We train a set multilevel manifolds (MLMs) using a database of visual hull reconstructions from spine x-ray images. We link to the constructed multilevel manifold, a set of parameters helping to generate the intervertebral articulations and shape parameters. This work significantly extends a preliminary version of the method, 13 where we demonstrated the potential of using multiple manifold instances to instantiate new models. However, this latter method was trained on preoperative models and used a charting function based on a mixture of Gaussian kernels, which was shown to be sensitive to outliers. Our main contribution is the development of an x-ray image spine reconstruction method, where the anatomy of the spine is occluded by instrumentation (screws and rods) as shown in Fig. 1. The method is based the development of an on-the-fly mapping process where local charts are generated whenever it is required on the manifold, estimating the tangent space in the vicinity of a data point and maximizing the accuracy for the nonlinear estimation. Our approach combines both building and learning the manifold structure within a neighborhood of this graph, in order to effectively infer new models.

MATERIALS AND METHODS
The proposed two-step procedure first learns an ensemble of multilevel manifolds from a training set of reconstructed spines obtained from a scoliotic population in order to capture the statistical variability of the geometrical shapes (Fig. 2). This ensemble includes multiple binary trees connected between each other to locate the nearest data points for any given new query point. The visual hull reconstruction of the global spine shape is used as the search descriptor and is obtained from a set of calibrated biplanar x-rays, using a silhouette extraction algorithm of the anterior portion of the spine extracted from both images. 14 In the testing phase, given a set of unseen biplanar x-rays, the visual hull is reconstructed to navigate toward the best match in the trained manifold ensemble. An image-based optimization procedure based on pairwise Markov random fields (MRFs) is finally applied to refine the 3D postoperative model. F. 1. Sample postoperative x-ray images from scoliotic patients following spinal surgery. Instrumentation includes metallic rods anchored on the vertebrae with pedicle screws attached through the spinous processes.
F. 2. Schematic representation of the 3D reconstruction process for postoperative spines. In the training phase, a multilevel manifold ensemble is learned from a dataset of 3D visual hulls partitioning the data. The model enables both within and in-between moves between tree nodes in the manifold datasets. The testing phase takes as an input a 3D visual hull, locates it within the manifold to perform nonlinear charting, and finally applies an image-based optimization procedure to generate the output 3D spine model.

2.A. Modeling the spine visual hull
In this first section, we present a method which has been previously described in fields such as in robotics, which takes advantages of the shape-from-silhouette principle in order to reconstruct the visual hull of the spinal column from posterior-anterior (PA) and lateral (LAT) radiographs. Toward this end, the segmented silhouette of the anterior section of the spinal anatomy obtained from a pair of calibrated x-rays is used for this purpose as it offers a global shape descriptor of the entire spine as shown in Fig. 3. The visual hull is represented by the pair of silhouettes s i with i = {PA,LAT}, and on the camera viewing parameters Π i (rotation, translation, and principal distances) such as Π i : ℜ 3 → ℜ 2 for the 3D to 2D projection, which is determined here using a conical or self-calibration technique. 14 Once the silhouettes are obtained from the images using a previously developed segmentation algorithm, 14 the reconstruction is based on the concept of visual hull. In other words, a visual hull of the 3D spine (denoted as x) based on a given projection plane Ω, represented as V H(x,Ω), can be defined as a 3D volume where for a given landmark l in V H(x,Ω) and for a given x-ray source V i in Ω, the projective vector from V i to l includes a point of x. 15 Therefore, we can assume that a visual hull object includes all points in 3D whose projections lie within the set of segmented shapes on the biplanar images from a radiographic system.
In this particular application, a set of corresponding segmented silhouettes as viewed in the projected x-rays are used to define the maximal object which represents the baseline spine model. The vertebral silhouettes are projected back into the 3D space by using a set of conical visual hulls in the projective model, with the projections coinciding with the spine outlines on Ω such that Π i (x) = s i . By determining the intersection of the conical projection obtained from the biplanar views, we can generate the estimated shape of the entire spinal column based on the calibrated parameters Π i . F. 3. Concept of the spine's shape 3D visual hull using the silhouette contours s PA and s LAT segmented on a pair of radiographic planes. The visual hull x is determined using the intersection of the conical projections originating from viewpoints V, yielding a high-level representation of the scoliotic spine.

2.B. Building the multilevel manifold ensemble
The training phase of the proposed method consists of learning the shape variations from a collection of visual hull reconstructions and create partitioning trees that randomly covers the training set of spine models. This concept follows previous approaches used to estimate human poses based on a multiresolution manifold. 12 Each partitioning tree is connected to all the other trees and can be defined as representing the ambient space using search grids that can adapt to the underlying data distribution. In this application, the multilevel manifold ensemble was trained from 543 reconstructions from postoperative scoliotic patients treated for various types of deformation. The postoperative reconstructions were generated using a semisupervised software yielding an initial parametric model. 4 Each model is then refined by an experienced user to adjust for anatomical landmark positions. For each spine in the dataset, both a visual hull reconstruction and a high-resolution geometrical model were obtained from calibrated biplanar x-rays. The models include 12 thoracic and 5 lumbar vertebrae, each represented by 10 landmarks (4 pedicle extremities, 2 endplate center points, and 4 anterior/posterior points on the superior and inferior endplates), yielding a total 170 landmarks per model.
One can define an ensemble of decision trees t i ∈ T that separate the ambient dataset R D of visual hulls. Given a dataset X = {x i }, x i ∈ R D of postoperative spine visual hulls of dimension D, the procedure trains a collection of decision trees composed of a set O of nodes. Here, we make the assumption that each 3D spine model lies on low-dimensional space M, of dimension d. Each node of a tree is defined by a separating hyperplane Φ i which is obtained from the higherdimensional space. Hyperplanes Φ i are represented by their unit vector normal n j defined in ambient space and a threshold value τ j which represents the average of two mapped sample points, which exhibits the largest Euclidean distance within the embedded space. The dataset of visual hulls is separated into two groups, namely, as X L j and X R j , which are determined based on a separating function defined as which avoids mapping data samples away from the manifolds. The set X L j represents points x ∈ X where the split function yields 0, while X R j contains samples where the split function yields 1. In order to determine the node parameters Φ i from the training data, a subset S j of visual hulls were randomly sampled from the general distribution of X j . From the subset S j , a single visual hull x k is selected and the data point further away (highest Euclidean distance) from x k within S j was determined. The normal n j of the high-dimensional plane is defined as the unit vector from the selected visual hull x k to the maximal point x max such that Once the partitioning trees are generated from the training data, the procedure learns the graphs from the nodes O and the set of directed edges, E. Here, we denote nodes as o t i ∈ O with t representing a given partitioning tree, while the node index given by i. An edge set is composed of all parent-child tree edges E t along with edges linking data structures E s, t . A collection of edges can be described as all the edges E t of trees t included in T , as well as the edges E s, t in the connecting trees s,t included in T . The rationale for the connected trees method is that two nodes o s i and o t j in separate trees can in fact be connected between each other if they show similar characteristics with their overlapping intersected region as shown in Fig. 4. Because the calculation of intersecting regions may be too difficult within the ambient space, the alternative approach seeks to estimate these regions if samples extracted from their neighboring are in fact shared between both trees.

2.C. Optimization on the manifold
Once the multilevel manifolds are trained, our objective, given a query case x q , is to optimize an energy function f which is determined from data points embedded on the manifold M. The first step seeks to obtain an approximate solution by navigating the entire multilevel manifolds, both through and across the multiple manifold instances (within and between trees). To facilitate the search process based on an optimized k-d search, the trees are aligned from the principal components in order to reduce the effect of variable rotation within the datasets. Once the desired leaf node has been reached, the local neighborhood of the node within the same tree is queried in order to generate a nonlinear chart on M.

2.C.1. Generating the approximate solution
For the minimization of f based on the search within the multilevel manifold framework, a collection of starting points F. 4. Illustration of connecting tree nodes that possess overlapping regions. Two separate trees describe different data partition trees that separate the overall ambient domain. Here, a multilevel manifold structure links various tree nodes where their respective regions of the dataset intersect. In this example, this is the case of the nodes included in both light and darker trees. Connecting tree nodes allow horizontal moves between trees. Data samples of neighboring data points are used to estimate intersections. are generated within the manifold space. These points offer a set of potential candidates that will fall near the desired solution. The goal of the procedure is not to generate a unique solution, but instead a collection of nodes that serve as a basis for the charting process that will create a nonlinear mapping function. Because the search process with the multilevel trees is rather complex, we employ a simplified approach to obtain the closest points from a query input. Prior to the tree search, the strategy will apply a PCA on the set of binary trees t i ∈ T in order to align them based on their respective principal components. The datasets are divided in the tree based on hyperplanes crossing the axis n j defined earlier, by first translating the data so that their centroids coincide at the origin. Then, we Here, the principal axis defined with the training set represents the eigenvectors of U. Consequently, one can therefore use the eigenvalues of U to represent the principal moments which are used to rotate the trees t i ∈ T .
For any new query point x q , a descent down the trees brings to a solitary node of the tree. A data instance linked with the node o t j corresponds to a first candidate in the search of the closest samples. Since this point may not be the closest neighbor, the process incorporates a search with horizontal and vertical moves, where additional nodes are visited to identify the best candidate points. Here, we adopt a priority search method 16 where an ordered sequence of tree nodes are searched based on their respective distance to the input data vector x q . By using a priority tree in order to reorganize the nodes, this can be achieved efficiently. In this implementation, tree nodes are associated with a specific number proportional to the Euclidean distance to x q . The process ends as soon as all tree nodes within a predefined range have been visited.
After the minimization procedure, the algorithm is able to generate a collection of pointsx t j representing the mean of a group of nodes x t j . This data serves as the basis to locally chart the manifold M around the query point x q . This ability represents a clear advantage to alternative technique based on graph-building, as the mapping function are determined on the fly over the training data.

2.C.2. Manifold mapping from kernel regression
During the navigation procedure on the manifold M, the nonlinear neighborhood is defined by leaf averagesx t j around a query point. Using a random walk including horizontal displacement (transitions between search trees), the seed points are propagated to their neighboring points. The search proceeds as long as b nodes are visited. The dimension of the neighborhood based on the value of b determines the extent of the local chart dimension and can vary based on the specific application and the optimization procedure that is employed.
In order to estimate the mapping from a local chart near x on the manifold M at y ∈ R d to the ambient domain in Dspace, we use the tree nodes obtained by the search process as samples to train the inverse mapping where the tangent plane serves as a basis to estimate the d-dimensional space and determine the nonlinear neighborhood charting. A preliminary version of this work used a Gaussian mixture models to perform the inverse mapping but was shown to be more sensitive to outliers. 13 Here, we propose to use a kernelized extension of multioutput ridge regression initially proposed for out-of-sample problems, 17 to allow for the mapping from y, representing the low-dimensional coordinates, to the ambient point x such that R d → R D . We use a Mercer kernel, which is a symmetric positive definite real-valued function inducing a map ψ that nonlinearly projects the manifold space M to a Hilbert feature space H, such that ⟨ψ( y i ),ψ( y j )⟩ = Q( y i , y j ). We denote here ⟨·,·⟩ the inner product in the space H.
We define Ψ = (ψ( y 1 ),ψ( y 2 ),...,ψ( y n )) as a low-dimensional vector defined in Hilbert space H, with W = [W 1 ,W 2 ,..., W d ] a vector of weight elements in H. The kernel ridge regression defined in feature space H can be defined as which produces the ambient vector from the projection of the low-dimensional coordinates and where ϵ is the vector of random residuals. Since the inner product ⟨Ψ,W⟩, defined in Hilbert space, is of infinite dimensions, this thereby reduces the matrix product to finite dimensions. The components W can be determined by the following minimization: where ∥ · ∥ F represents the Frobenius norm of a matrix and λ is a regularization parameter which finds the best compromise between the project error and the penalty term. The solution of problem (4) is determined explicitly by setting the first derivative to zero. By derivation, the predicted representation x for a low-dimensional point y can therefore by found as where I represents the identify matrix, κ( y) = ⟨Ψ,ψ( y)⟩ = (Q( y 1 , y),Q( y 2 , y),...,Q( y d , y)) represents the kernel function used for the inverse mapping, and K is the Gram matrix such that K i j = Q( y i , y j ) based on the kernel mapping function Q. This constrains the regression to be valid for all the neighborhood of the manifold chart around y i as it preserves locality in x i . This nonlinear space parameterization is only valid for a predefined tangent space. For any new space domain, a new ridge regression must be obtained within the neighborhood of the updated solution.

2.D. 3D reconstruction of postoperative patients
In order to reconstruct 3D models of the spine from biplanar images where the anatomical landmarks are partially occludes with metallic instrumentation, we describe the process as a navigation on top of the manifold of visual hulls from previously reconstructed instrumented spines. The objective function formulated here looks to obtain the best match between the input visual hull of the postoperative spine and observations in the manifold. The smoothed visual hull shapes, represented as x, are used to navigate in the manifolds based on the query function f based on Ref. 12 in order to find similar instances to x, such that (6) which describes the similarity between the input x and points y on the manifold. The kernel mapping k(·) is able to limit the impact of potential outliers, and I(·) indicates whether the input vector points x i fall within a certain range ϵ to the manifold points y i . The function f incorporates a penalty term β for samples that lie far away from the input reconstruction x or that fall inside the visual hull, i.e., the instrumented visual hull must lie as close as possible to the input shape. Once approximate models are generated from the manifold, the regression process in Eq. (5) defined over a local region enables to refine the geometrical model of the postoperative spine.
Given a reconstructed model from the hierarchical manifolds, the reconstruction process reverts back to the highdimensional domain using the mapping function φ( y i ) in Eq. (5) to generate the complete geometric model. Final refinements are introduced to the geometrical model by improving the fit of the retro-projected model onto the biplanar xrays, using a 3D-2D optimization paradigm. 18 We build the spine geometry from the generating parameters linked to the manifold neighborhood determined previously, depicting the spine by an array of local intervertebral rigid registrations A = [T 1 ,T 2 ,...,T m ], with m the number of intervertebral transformations (m = 16) and T i = {R,t} a rigid transform (rotation and translation) between consecutive vertebrae. We use a higherorder graphical model to adjust vector A, which was obtained by the reverse manifold mapping process. Here, the graphical model is described by a MRF which includes an image similarity term g = {g i (·)}, defined for all nodes in the MRF, a pairwise potential p = {p i j (·,·)} acting as a regularization term to constraint displacements between nodes and a higher-order term c = {h c (·)}, which denotes the higher order potentials, representing the displacement of nodes within a clique c. The data term describes geodesic active contours (GACs) from the 3D vertebra models projected into the gradient of the radiographic images, enhancing the vertebral cortical surface of the vertebral bodies. This was shown to effectively fit a volumetric data onto projected radiographic images. 7 The pairwise potential serves to limit important displacements between consecutive vertebrae (e.g., T 4 and T 5 ). This prior term helps to limit the displacements between pairs of local vertebral coordinates during minimization, by generating a pairwise potential p i j = {0,1} such as with ℜ the transformation field of a node (ℜ ∈ t,R), ⃗ d i and ⃗ d j are displacements vectors associated to the transformations, while ϵ 1 ,ϵ 2 are empirically set values which were determined for specific vertebrae. The first condition restrains the translation within the same approximate direction. The other constraint imposes similar ranges of translation and rotation for a specific vertebral level. For example, a high deformation in translation usually involves high degree of rotation in scoliotic spines. This parametric modeling introduces coherency when minimizing the energy term of the graphical model. Finally, the higher-order term h c helps to regularize the displacement between triplets of vertebrae, specifically the apical (1) and two end vertebrae (2) for an anatomical region (thoracic, lumbar, proximal-thoracic). Additional details for this term can be found in previous work. 19 The MRF cost function to minimize is defined as Here, labels u i ,u j ∈ L are associated to transformations T i ,T j ∈ A, with T i modified by adding a translation or rotation as defined in the label space L, with u i ∈ L. On the other hand, u c ∈ L is linked to three transformations T i ,T j ,T k ∈ A, with {i, j,k} representing the apical and two end vertebrae in A. The energy term, defined in terms of the associated labels, is minimized using a discrete optimization procedure based on primal-dual decomposition and uses a nonsubmodular term to solve the above-defined MRF. 20

EXPERIMENTAL RESULTS
The reconstruction approach was evaluated based on a series of numerical, synthetic, and clinical experiments in order to quantify the precision of the manifold-based statistical model pipeline, as well as the overall 3D reconstruction accuracy of the method. A simulated 3D spine model was first produced, with varying levels of noise and occluded vertebral levels added to the model, allowing to set the appropriate parameter settings of the algorithm and study the robustness of the 3D reconstruction toward different initial conditions. The new reconstruction technique was also tested on synthetic spine models (reference) and real patient data, from both preoperative and postoperative images. Experiments were performed on a i7 Intel PC workstation, with 32.0 GB of memory.

3.A. Simulations on 3D computerized models
We first tested the multilevel manifold structure for the 3D reconstruction of spine models with or without instrumentation by testing the framework under controlled settings. A specific computerized spine model with a right thoracic deformity was created in a predefined coordinate system. In order to assess the general robustness and precision of the 3D models, simulations were undertaken by testing various conditions with regards to Gaussian noise and occlusions from x-rays. The computerized model was projected on a set of images using a typical x-ray setup in order to generate ideal conditions. For these experiments, an exhaustive search of parameter settings was performed to determine the optimal set of values. The node neighborhood size was b = 6, the Gaussian kernel size for the nonlinear regression model was σ = 1.5, and the penalty value β for points falling outside the visual hull was β = 0.8. As for the threshold values ϵ 1 ,ϵ 2 for the pairwise potentials, they were set to the same values as presented in previous work. 19 The overall robustness of the multilevel manifold ensemble model was first evaluated toward noise added to the input visual hull reconstruction. We tested the algorithm's ability to retrieve the closest match from the training dataset by comparing the mean 3D distance of the model to its ground-truth based on the accuracy of the 3D visual reconstruction. In this test, increasing levels of noise (11), ranging from 0 to ±5 mm, were added to the input visual hull meshes. In sum, 100 trials were performed for each noise level. Figure 5(a) presents the evolution of the 3D RMS distance of landmarks positions and surface mesh differences with respect to the noise level and demonstrates the capacity of the multilevel manifolds to locate the accurate model from the training dataset. In order to generate realistic configurations to surgical cases, spinal instrumentation masks were generated and added to an increasing number of vertebral levels along the length of the silhouettes in order to simulate spinal instrumentation. Based on the results presented in Fig. 5(b), the algorithm seems to be relatively sensitive toward occlusion near the upper-thoracic vertebrae which are smaller in size, thus more affected by occlusion from instrumentation as they hinder the visibility of the vertebral bodies. On the other hand, 3D reconstruction errors remain stable as they are added toward the lumbar region. Results demonstrate the framework's robustness toward noisy input data and illustrates that adding regularization potentials allows to capture the overall distribution of anatomical landmarks.
Finally, we evaluated the convergence of the 3D reconstruction process per iteration using three different methods for building statistical models. Figure 5(c) presents the convergence results for this experiment. The MLM approach using nonlinear mapping during navigation, out-performing global PCA and a kernel LLE method 21 for the initialization of the model using the extensive database of models, as well as the later optimization stage with image support to refine the vertebral models on the x-rays.

3.B. Accuracy testing on synthetic models
The accuracy of the method was evaluated based on our previous methodology 14 using an in vitro reference spine model with three different configurations of scoliosis: normal, rightthoracic left-lumbar, and lumbar deformity. A specific setup offers a representation of a particular class of deformation which can be seen in the scoliotic population. In all, 19 spher-  ical markers of 2 mm in diameter were embedded in each vertebra of the synthetic model. These pellets were placed on the vertebral body and spinous processes. A 3D digitizer (MicroScribe, accuracy of 0.3-mm) was used to measure the precise 3D position of each pellet in the global coordinate system for registration and comparison purposes (Fig. 6). After measuring these coordinates, a rigid registration (global and local) was applied in order to transpose the digitized spine in the same coordinate system as the biplanar 3D reconstruction by using 5 markers placed on the spinous processes, transverse processes and laminae. Absolute position and orientation of the 3D reconstructions are not important since what must be measured is how accurately the general shape of the spine is reproduced.
The accuracy of the 3D landmark positions generated by the multilevel manifolds was obtained by measuring the pointto-point Euclidean distances of the vertebral body landmarks (n = 10), as well as the superior and inferior pedicle extremities (n = 4) between the 3D digitized model and the generated reconstruction. With steel pellets of known 3D coordinates incrusted in each vertebra of the spine, this enables us to quantify the reconstruction error.
The 3D mean distance between the models generated from the 3D reconstruction method and the measured ground-truth is presented in Table I. Overall, the error was of 2.0 ± 1.0 mm for endplate landmarks and 3.0 ± 1.5 mm for pedicle landmarks, yielding a global error of 2.3 ± 1.2 mm for all landmarks. As previously shown, 14 a possible source of error can be accounted from the transformation errors from the local to global coordinate systems, which was globally evaluated at 0.3 ± 0.15 mm. However this level of discrepancy is suitable to evaluate the general posture of the spine and allow for clinical diagnostic/evaluation applications. While the previous simulated experiments enable to test the method's tolerance to external factors, the accuracy level presented from these synthetic experiments compared favorably to previous 3D reconstruction methods. 22

3.C. Preoperative scoliotic spine reconstruction
We then evaluated the performance of the method on preoperative patients with scoliosis. Ten pairs of preoperative biplanar x-rays taken from patients with mild deformities (Cobb angle 15 • -55 • ) were used to quantify the 3D accuracy of the T I. Mean and RMS accuracy of the reconstruction framework obtained on three synthetic spines digitized in vitro with high precision. The geometrical coordinate results of the markers are presented at each vertebral level (T1 to L5) for each anatomical group (endplates and pedicles).

Endplate markers (mm)
Pedicle markers (mm) Global error (mm) generated models. This enables to assess the errors under normal settings, without any occlusion from instrumentation. For each case, differences between the proposed method and landmarks identified manually by a radiology expert were computed. Sample results are shown in Fig. 7. These promising results show the projected 3D vertebral shapes obtained from the reconstruction technique integrating multilevel manifold ensemble following the apparent structural edges of the vertebra's cortical surface on the pair of x-ray images. Table II presents the quantitative results, demonstrating performances similar to state-of-the-art, using a single manifold approach. 7 With d = 10, the mean difference in 2D landmark locations for all ten patients was of 1.96 ± 0.65 mm for thoracic vertebrae and 1.25 ± 0.39 mm for lumbar vertebrae, while the 3D RMS difference was 2.20±0.74 mm for thoracic vertebrae and 1.32 ± 0.47 mm for lumbar vertebrae. While the improvement of the method in terms of 3D reconstruction errors compared to the previous approach based on locally linear embedding was not significant using a paired student t-test analysis (p < 0.05), there is a consistent reduction in the 3D RMS difference compared to manual annotations.

3.D. Postoperative spine reconstruction
Finally, a testing set of thirty (n = 30) pairs of instrumented biplanar x-rays were used for postoperative reconstruction, where 3D models were generated by learning a manifold of visual hulls of spine models. Sample results are shown in Fig. 8. In Table III 35 mm), based on a paired student t-test analysis (p < 0.05), as well as an improvement to a preliminary version of this work. 13 To evaluate robustness of the method toward different instrumentation strategies, we measured the reconstruction accuracy over five different types of surgical strategies, based on instrumented levels (thoracic and lumbar). These include upper-thoracic (n = 7), mid-thoracic (n = 9), full thoracic (n = 5) (upper and mid sections), lumbar (n = 4), and thoraciclumbar (n = 5) regions. Since various surgical strategies will cover different anatomical regions in the spine, this experiment can demonstrate the ability to handle occlusion at different levels of the spine. Figure 9 shows the encouraging 3D error results with different instrumentation levels. It also illustrates the benefit of integrating horizontal moves in the ensemble search when different anatomical regions are occluded by metal rods, allowing an efficient navigation toward the best match in the dataset.

DISCUSSION AND CONCLUSION
In this paper, we proposed an unsupervised method to perform the 3D reconstruction of a spine's geometry when surgical instrumentation is visible on biplanar x-ray images. Our approach is based on an ensemble of multilevel manifolds, which enable an efficient navigation on a low-dimensional domain and infer the closest match within a training set of previously reconstructed spine models from a population of scoliotic patients. An image-based optimization procedure based on higher-order Markov random fields is then used to refine the spine model based on the fit of the vertebral shapes with the x-ray images, paired with regularization potentials to constrain the displacement between vertebrae. Results show that this model allows reconstruction accuracies similar to gold-standard. Because the 3D reconstructions are obtained with little user interaction, the proposed approach could be transposed to clinical practice and used in the context of multicenter evaluations of surgical practices and longitudinal evaluation of scoliosis progression.
One of the challenges in reconstructing the spine in 3D from postoperative x-ray images is poor visibility due to the occlusion caused by the spinal instrumentation, such as the metallic rods and screws, as well as the superposition of multiple anatomical structures. These problems are particularly apparent in the lateral view near the rib-cage area. Typically, an operator or radiologist has to rely on past experience and knowledge of the anatomy in order to estimate the landmark locations on the vertebrae and visually infer landmark positions even when rods and screws are obstructing the appearance of the anatomy. The proposed framework based on the navigation of multilevel manifolds increases the overall reliability of the reconstruction approach by including statistical knowledge of the spine appearance variability. This allows to offer a suitable platform for the reconstruction of spine models for postoperative evaluations. The experimental evaluation shown here was crucial to portray the reliability of the multilevel manifold ensemble to accurately predict the position and rotation for constelations of vertebrae with increasing levels of complexity that are usually observed in clinical practice. The method yields 3D reconstructions that were sufficiently close to models obtained by an expert operator. Unfortunately, supervised approaches are cumbersome and sensitive to errors, and cannot guaranty 100% accuracy. Thus, differences present in this experiment may originate from identification errors generated by the expert user.
Radio-opacity from the metallic instruments and varying levels of x-ray doses can affect the appearance of the vertebral anatomy, generating suboptimal visual hull reconstructions that can ultimately affect the navigation within the MLM data structure. Furthermore, input reconstruction errors within the training spinal geometries can affect the training phase of the manifold ensemble and ultimately change the behavior of the predictive model. Still, our experiments show the method reconstructs the spinal shape with high accuracy even with blurred or occluded contours caused by the instrumentation artifacts. This can be attributed to the inclusion of regularization constraints which prohibits inconsistent modification to the initial geometry. An advantage of automated reconstruction methods lies in the increased repeatability for identifying anatomical landmarks which reflects the statistical distribution from prior reconstructions. This helps to reduce the intrinsic variability induced by manual localization by effectively capturing the global (overall spine shape appearance) as well local (vertebral representation) within a study group, thereby offering higher tolerance to patient variability. The model can also handle out-of-plane 3D rotation such as vertebral axial rotation, which in the case of AIS, is frequently observed by clinicians and where manual techniques exhibit a number of limitations. 23 The potential impact of this 3D reconstruction technique was shown by treating a group of preoperative and postoperative scoliotic patients. Results confirm that using a multilevel manifold structure to search for the ideal candidate used for the reconstruction process, coupled with an optimization procedure based on local shape refinement, enables the precise reconstruction of the spine geometry. In fact, the visual hull representation of the spine can serve as a reliable highlevel primitive to automatically navigate through the manifold structure. Indeed, this general descriptor offers rich features that can be easily exploited by the MLM structure in order to partition the data in multiple tree instances. The proposed automatic technique allows generating accurate 3D spinal shapes compared to other dimensionality reduction techniques based on linear projections (PCA) or standard nonlinear techniques (LLE) that are less adapted to handle highly variable representations and out-of-sample problems. By comparing the 2D and 3D RMS reconstruction errors with the proposed method and the statistical modeling techniques method, we illustrated that the convergence of the manifold structure obtains the lowest levels of error, in addition of reaching the global optimum with less iterations. Indeed, the proposed approach includes an efficient manifold navigation method that integrated a nonlinear mapping process which infers new models from local neighborhoods within the different manifolds. Unlike traditional manifold learning techniques, the multiple manifolds include not only top-down links to allow to locate the closest model from the training set, but horizontal connections between the multiple manifolds, which enables a more refined search of the model. Furthermore, the image-based optimization procedure includes parametric constraints that improve regularization. While the accuracy of the method is promising, errors may be propagated from the visual hull extraction due to potential patient motion during the radiographic acquisition, ultimately affecting the algorithm's convergence since the pairwise correspondence between the images is affected. Still, it allows a reconstruction process of postoperative scoliotic spines, previously very difficult to achieve using methods that require a high level of manual intervention and outperforms a previous version of the work, 13 by training the manifolds with datasets of instrumented cases thus improving the search process, as well as proposing a more efficient charting process based on ridge regression.
Future work will improve this approach by including morphologic parameters, as well as neighborhood relationships between the structures to increase the reliability and repeatability of the reconstruction approach. These parameters can be modeled in a top-down representation of the spine anatomy and representing the intervertebral constellation within an articulated model. 24 Currently, the training is performed on an extensive database of preoperative spine models, due to the limited amount of postoperative reconstructions available. We plan to improve our approach by using a more adapted training database to capture the variability within the population. The approach can be used for other shape modeling problems, with partial view challenges, for example, with the distal femur in orthoplasty procedures, and adapted to different applications in medical imaging for 3D shape reconstruction problems (fluoroscopic tracking), or in computer vision for scene inference.