Reconstruction of three-dimensional lumbar vertebrae from biplanar x-rays

Three-dimensional models that are reconstructed from medical images have significantly aided the process of clinical diagnosis and surgical planning because they can visualise the internal musculo-skeletal structures and tissues with high geometric accuracy without open wounds. Common sources include computed tomography (CT), magnetic resonance imaging (MRI), and ultrasound. CT scan reconstructs both bony and soft tissues, but the constant radiographic scan around the body for producing tomograms would subject the patient to high doses of radiation. On the other hand, MRI, which is suitable for visualising soft tissues, is considered a safer option because it utilises an oscillating magnetic field to generate tomograms, but patients with pacemakers and metallic implants may not undergo an MRI exam safely. Ultrasound utilises high-frequency sound wave to produce tomograms (Kak and Slaney 2001), but the reconstructed model can be noisy.

Orthopaedic surgical navigation systems rely on CT-reconstructed spine models for planning and guiding. These systems, coupled with optical tracking systems and robots, can deliver percutaneous operations with minimum risks (Hyun et al 2017, Khan et al 2019, Mao et al 2020). These systems require engineers to setup, and pre-operative CT scans may expose excess radiation to the patient. There are fluoroscopic-based guiding and navigation systems that takes intra-operative radiographs in anterior-posterior (AP) and lateral (Lat.) views to plan and display tool paths, but extrapolating 2D information to 3D can be difficult (Foley et al 2001). It would be beneficial to reconstruct vertebrae models from the radiographs, as it reduces the of exposed radiation dosage compared to CT reconstruction, and visual reference on the vertebrae's relative positions can be obtained.

The earliest work in visualising spinal deformities from biplanar x-rays employed analytical photogrammetry to establish calibration between two biplanar x-rays using collinearity equations, with lead targets attached to the x-ray cassette holder serving as the calibration object (Kratky 1975). The effect of the chosen vertebral landmarks, viewed from two stereographic planes, was studied on the 3D reconstruction accuracy, and had established six stereo-corresponding landmarks that produced the least error for a vertical stereographic setup (André et al 1994); however, the outcome did not include local geometric details.

To account for this, the non-stereo-corresponding point (NSCP) method was developed by Mitton et al (1994), in which besides the six stereo-corresponding points, points that are only visible on a single radiograph were also identified to describe the local features. It yielded a mean error of 1.1 mm for cadaveric lumbar vertebrae, using 19 NSCPs on average (Mitulescu et al 2001), and was further validated on scoliotic patients, yielding a mean error of 1.5 mm (Mitulescu et al 2002).

Semi-automated methods using statistical models were introduced later to address the issue of extensive manual operation, and they achieved a high level of reconstruction accuracy comparable to CT models thanks to the a priori statistical shape models (Kadoury et al 2009). These methods required little user input (Kadoury et al 2009, Humbert et al 2009, Moura et al 2011), and rely on the information in the images such as vertebral contours (Kadoury et al 2009, Zheng et al 2010, Zhang et al 2013) or morphological features (Kadoury et al 2009, Stern et al 2013); some methods could reconstruct the model using only a lateral x-ray (Zheng et al 2010). These methods rely on a large statistical database of non-scoliotic and scoliotic spine specimens to parameterise geometrical characteristics.

The smoothness of the reconstructed model from template-based methods, such as the NSCP method, depends of the number of points; statistical methods can produce a more realistic reconstruction, but they can take longer to compute, and the limited information from the radiograph may hinder reconstruction accuracy (Fang et al 2019).

In some orthopaedic surgeries, multiple radiographs are taken by the surgeon to locate the inserted instrument and its relation to the target inside the body. Surgical experience and knowledge of spinal anatomy is crucial; less experienced surgeons can take longer to adjust the inserted needle to the target, and the extra fluoroscopic shots taken for confirmation would increase radiation exposure to the patient and surgical staff. Commercial percutaneous image-guided navigation systems are available, but they are more expensive, and CT-reconstructed models that are essential for the systems may be unavailable due to insufficient time and expense. A method that generates patient-specific lumbar vertebrae models during operation for the surgeon's reference while minimising radiation exposure is required.

We therefore proposed a lumbar vertebrae reconstruction method from biplanar intra-operative x-rays to generate patient-specific vertebrae models. The proposed method utilised non-stereo-corresponding landmarks, defined on AP and Lat. X-rays, to deform a set of lumbar vertebrae template models.

As depicted in figure 1, the use of the two x-rays is similar to a stereo camera system, where the coordinates of the landmarks on the vertebrae can be obtained and to be used for deformation. Figure 2 outlines the overall procedure. AP and Lat. X-rays of the region of interest, taken from a mobile C-arm, were transferred to our developed programme to be calibrated and corrected, and then non-stereo-corresponding landmarks were identified by selecting the vertebral landmarks on the images; the landmarks were then retro-projected to the 3D coordinate system. A set of lumbar vertebrae template models, as well as the corresponding landmarks on them, were imported, and their orientation was estimated. The template models that corresponded to the identified vertebrae were deformed by matching the template landmarks to their corresponding user-identified landmarks using free-form deformation driven by an optimisation problem, yielding a set of deformed patient-specific lumbar vertebrae models. The programme could identify individual landmarks based on their relative locations; the operator would be notified on any mis-matched landmarks from the AP and Lat. X-rays to prevent the programme from producing spurious models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. Camera calibration and distortion correction

The principle of image formation of a C-arm XRII is like that of a conventional pin-hole camera model, which establishes the homography that sets up the linear transformation between 2D image coordinates ${\left[uv1\right]}^{{\rm{T}}}.$ to 3D space coordinates ${\left[XYZ1\right]}^{{\rm{T}}}.:$

$$\begin{eqnarray}&&\begin{array}{c}s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=A\left[{\boldsymbol{R}}\left|{\boldsymbol{t}}\right.\right]\left[\begin{array}{c}\begin{array}{c}X\\ Y\end{array}\\ \begin{array}{c}Z\\ 1\end{array}\end{array}\right]\end{array}\end{eqnarray} \tag{ 1 }$$

where A is the camera intrinsic parameters, and $\left[{\boldsymbol{R}}\left|{\boldsymbol{t}}\right.\right]$ is the camera extrinsic parameters; s is the size of a pixel on the C-arm XRII, which was set to 1.

To remove distortion effects on the x-rays, they were fitted to a fourth-order polynomial (Fahrig et al 1997) to compute the global distortion coefficients; higher-order polynomial may result in over-fitting. Let the distorted coordinates on the image be (U_n , V_n ), and their corresponding non-distorted coordinates be (u_n , v_n ), their relationship can be expressed as:

$$\begin{eqnarray}\begin{array}{l}\left[\begin{array}{cc}{u}_{1} & {u}_{1}\\ {v}_{1} & {v}_{1}\end{array}\ldots \begin{array}{c}{u}_{n}\\ {u}_{n}\end{array}\right]\left[\begin{array}{cc}{a}_{0} & {a}_{1}\\ {b}_{0} & {b}_{1}\end{array}\ldots \begin{array}{c}{a}_{14}\\ {b}_{14}\end{array}\right]\\ =\,\left[\begin{array}{ccccccccccccccc}1 & {U}_{1} & {V}_{1} & {U}_{1}^{2} & {U}_{1}{V}_{1} & {V}_{1}^{2} & {U}_{1}^{3} & {U}_{1}^{2}{V}_{1}^{2} & {U}_{1}^{2}{V}_{1} & {V}_{1}^{3} & {U}_{1}^{4} & {U}_{1}^{3}{V}_{1} & {U}_{1}^{2}{V}_{1}^{2} & {U}_{1}^{2}{V}_{1}^{3} & {V}_{1}^{4}\\ 1 & {U}_{2} & {V}_{2} & {U}_{2}^{2} & {U}_{2}{V}_{2} & {V}_{1}^{2} & {U}_{2}^{3} & {U}_{2}^{2}{V}_{1}^{2} & {U}_{2}^{2}{V}_{2} & {V}_{2}^{3} & {U}_{2}^{4} & {U}_{2}^{3}{V}_{1} & {U}_{2}^{2}{V}_{2}^{2} & {U}_{1}^{2}{V}_{2}^{3} & {V}_{2}^{4}\\ \, & \, & \, & \, & \, & \, & \, & \, & \vdots & \, & \, & \, & \, & \, & \,\\ 1 & {U}_{n} & {V}_{n} & {U}_{n}^{2} & {U}_{n}{V}_{n} & {V}_{n}^{2} & {U}_{n}^{3} & {U}_{n}^{2}{V}_{1}^{2} & {U}_{n}^{2}{V}_{n} & {V}_{n}^{3} & {U}_{n}^{4} & {U}_{n}^{3}{V}_{1} & {U}_{n}^{2}{V}_{n}^{2} & {U}_{1}^{2}{V}_{n}^{3} & {V}_{n}^{4}\end{array}\right],\end{array}\end{eqnarray} \tag{ 2 }$$

in which the distortion coefficients, ${a}_{i}$ and ${b}_{i},$ can be solved using the least-squares method.

2.2. Landmarks identification

We identified non-stereo-corresponding landmarks on the x-rays ('target landmarks') to be projected to the 3D coordinate system; their corresponding 3D landmarks ('template landmarks') on the template vertebrae models, as well as the models' orientation in space, were also determined, before commencing model deformation.

The target landmarks could only be observed from a single projection, and from the anatomical perspective, they correspond to the anatomies the surgeon would identify and picture in 3D from the biplanar x-rays during orthopaedic operations. A total of 16 landmarks per lumbar vertebra were chosen to be the basis of deformation, with 10 landmarks on AP x-ray, and 6 landmarks on Lat. X-ray (figure 3):

• AP landmarks include the corners of the vertebral body, the upper and lower tips of each transverse process, and the upper and lower tips of the spinal process;
• Lat. landmarks include the corners of the vertebral body and the upper and lower tips of the spinal process.

Zoom In Zoom Out Reset image size

The spinal process landmarks on AP x-ray were for determining the roll angle; they were not directly involved in model deformation.

The template landmarks were pre-determined with a set of template vertebrae, and they can be used repeatedly for the same set of template vertebrae. On the template vertebrae, the transverse process and the vertebral body posterior landmarks were located on its vertices; they could be easily found on the corner vertices of the anatomies' template.

The other template landmarks were located on designated planes that intersected the template: the template's optimal symmetry plane (OSP) (Wong et al 2005) and mid-body plane. The OSP is a plane from which the most voxel pairs across it can be found, and it can be used to find the vertebral body anterior landmarks and the spinal process landmarks. The mid-body plane lies approximately across the widest section of the template vertebral body; the procedure to determine the mid-body plane is shown in figure 4. The vertebral body AP landmarks can be found on it.

Zoom In Zoom Out Reset image size

The approximate locations of landmarks were found before determining their true locations. The triangular facets that intersect either the OSP or the mid-body plane were first identified. Their vertices were projected to $yz$ and $xz$ planes to obtain their respective centroids, and then the projected vertices were divided into four quadrants. From the facets that intersected the OSP, the two furthest vertices from their centroid in the second and third quadrants were the approximate landmarks of the vertebral body anterior landmarks; for the mid-body plane, the furthest vertices in all quadrants were deemed as the approximate landmarks for the body's AP-view landmarks. Then, on each of the facets that the approximate landmark belongs to, we calculated the landmark's shortest path to the designated planes. In the end, the overall shortest path was used to determine the true landmark location, which was on the intersection between the path and the intersecting line between its corresponding facet and the designated plane.

The problem of finding the true landmark location from the approximated landmark vertex is stated below. Let a plane ${\rm{\Pi }}$ intersects a model triangle facet ${{\rm{\Pi }}}_{f}$ with a normal vector $\bar{n}$ at two vertices ${v}_{{\rm{i}}1}$ and ${v}_{{\rm{i}}2},$ and the goal is to find vertex ${v}_{{\rm{true}}}$ on ${\rm{\Pi }}$ that is the closest to the approximated landmark vertex ${v}_{{\rm{approx}}}$ and is on the facet plane ${{\rm{\Pi }}}_{f}.$ The shortest path $\overline{{l}_{p}}$ between ${v}_{{\rm{approx}}}$ and ${v}_{{\rm{true}}}$ would also lie on ${{\rm{\Pi }}}_{f}.$ Point ${v}_{{\rm{true}}}$ is located somewhere on the line $\mathop{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\limits^{\longleftrightarrow},$ and can be found by finding the intersection between the shortest path $\overline{{l}_{p}}.$ from ${v}_{{\rm{approx}}}$ to line $\mathop{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\limits^{\longleftrightarrow}.$ The shortest path from a point to a line is always perpendicular to the line; the shortest path $\overline{{l}_{p}}$ can be determined from the cross product between $\bar{n}$ and $\overline{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}:$

$$\begin{eqnarray}&&\begin{array}{c}\overline{{l}_{p}}={v}_{{\rm{approx}}}+\left(\bar{n}\times \overline{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\right)\cdot t\end{array}\end{eqnarray} \tag{ 3 }$$

Where $t$ is a scalar value. Point ${v}_{{\rm{true}}}$ can then be determined by finding the intersection between $\overline{{l}_{p}}$ and $\mathop{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\limits^{\longleftrightarrow}:$

$$\begin{eqnarray}&&\begin{array}{c}{v}_{{\rm{approx}}}+\left(\bar{n}\times \overline{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\right)\cdot t={v}_{i1}+\overline{{v}_{{\rm{i}}1}{v}_{{\rm{i}}2}}\cdot u,\end{array}\end{eqnarray} \tag{ 4 }$$

where $t$ and $u$ are scalars to be solved.

2.3. Rotation of template

2.4. Template vertebrae deformation

The free-form deformation (FFD) method (Sederberg and Parry 1986) was implemented to deform the vertebrae templates, which was driven by an optimisation problem. The geometry of the deformed model was defined by the trivariate Bernstein polynomial function.

$$\begin{eqnarray}&&\begin{array}{c}{p}_{{\rm{FFD}}}=\displaystyle \sum _{i=0}^{l}\displaystyle \sum _{j=0}^{m}\displaystyle \sum _{k=0}^{n}{B}_{i}\left({p}_{u}\right){B}_{j}\left({p}_{v}\right){B}_{k}\left({p}_{w}\right){v}_{ijk}\end{array}\end{eqnarray} \tag{ 5 }$$

where ${p}_{u},$ ${p}_{v}$ and ${p}_{w}$ specify the location of the point inside the deformation lattice; ${p}_{{\rm{FFD}}}$ is the location of a point after deformation; ${v}_{ijk}$ are the control vertices; and ${B}_{i}\left({p}_{u}\right),$ ${B}_{j}\left({p}_{v}\right)$ and ${B}_{k}\left({p}_{w}\right)$ are Bernstein polynomials in degrees $l,$ $m$ and $n,$ respectively, for example:

$$\begin{eqnarray}&&\begin{array}{c}{B}_{i}\left({p}_{u}\right)=\left(\begin{array}{c}l\\ i\end{array}\right){\left(1-{p}_{u}\right)}^{l-i}{p}_{u}^{i}\end{array}\end{eqnarray} \tag{ 6 }$$

where $\left(\begin{array}{c}l\\ i\end{array}\right)$ is the binomial coefficient. The deformation lattice was defined as each template vertebra's object-oriented bounding box (OOBB). By observing the distribution of the landmarks on the template model, the number of planes –$l,$ $m$ and $n$—in each direction were designated to $l=2,$ $m=6,$ and n = 2, resulting in 24 control vertices. This was so that we could control the movement of each landmark independently with designated control vertices, while keeping the number of control vertices as low as possible.

The goal of model deformation was to match each template landmark to its corresponding target in its view. That is, we tried to minimise the sum of distances $d$ between each displaced template landmark—${p}_{i}^{{\rm{AP}}}$ and ${p}_{i}^{{\rm{LAT}}}$—and its corresponding target—${t}_{i}^{{\rm{AP}}}$ and ${t}_{i}^{{\rm{LAT}}}$—respectively in its view. We defined an optimisation problem to approach this goal. There are 72 variables, ${\boldsymbol{x}}={x}_{1},{x}_{2},\ldots ,{x}_{72},$ to the problem, representing the movements of the 24 control vertices: ${x}_{1}$ to ${x}_{24}$ represent the movement in the local coordinate's $\vec{u}$ direction; ${x}_{25}$ to ${x}_{48}$ represent the movement in $\vec{v}$ direction; and ${x}_{49}$ to ${x}_{72}$ represent the movement in $\vec{w}$ direction. By expressing each template landmark (${p}_{i}^{{\rm{AP}}}$ and ${p}_{i}^{{\rm{LAT}}}$) in terms of the trivariate Bernstein polynomial, as in equation 4, the objective function can be formulated as follow:

$$\begin{eqnarray}&&\begin{array}{c}Min\,d=\displaystyle \sum _{i=1}^{{n}_{{\rm{AP}}}}{\left({\left({t}_{i}^{{\rm{AP}}}.x-\displaystyle \sum _{j=1}^{{n}_{{\rm{v}}}}\left({x}_{j}{B}_{ij}^{{\rm{AP}}}\right)\right)}^{2}+\left({t}_{i}^{{\rm{AP}}}.z-\displaystyle \sum _{j=1}^{{n}_{{\rm{v}}}}{\left({x}_{j+48}{B}_{ij}^{{\rm{AP}}}\right)}^{2}\right)\right)}^{\displaystyle \frac{1}{2}}\\ +\displaystyle \sum _{i=1}^{{n}_{{\rm{Lat}}}}{\left({\left({t}_{i}^{{\rm{Lat}}}.y-\displaystyle \sum _{j=1}^{{n}_{{\rm{v}}}}\left({x}_{j+24}{B}_{ij}^{{\rm{Lat}}}\right)\right)}^{2}+{\left({t}_{i}^{{\rm{Lat}}}.z-\displaystyle \sum _{j=1}^{{n}_{{\rm{v}}}}\left({x}_{j+48}{B}_{ij}^{{\rm{Lat}}}\right)\right)}^{2}\right)}^{\displaystyle \frac{1}{2}}\end{array}\end{eqnarray} \tag{ 7 }$$

Subjected to

$$\begin{eqnarray}&&\begin{array}{c}{x}_{k}={x}_{k}0,\,k=1,2,9,10,11,12,13,14,21,22,23,24\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{c}{x}_{k}={x}_{k+1},\,k=25,27,29,\ldots ,47\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\begin{array}{c}{x}_{k}={x}_{k+1},\,k=49,57,59,61,69,71\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{c}{{\boldsymbol{x}}}_{0}-50\leqslant x\leqslant {{\boldsymbol{x}}}_{0}+50\end{array}\end{eqnarray} \tag{ 11 }$$

where ${n}_{{\rm{AP}}}$ and ${n}_{{\rm{Lat}}}$ are the number of AP and Lat. landmarks, respectively; ${n}_{{\rm{v}}}$ is the number of control vertices; ${B}_{ij}^{{\rm{AP}}}$ and ${B}_{ij}^{{\rm{Lat}}}$ are the trivariate Bernstein polynomial functions of the ${i}^{th}$ template landmark in AP and Lat. view, ${p}_{i}^{{\rm{AP}}}$ and ${p}_{i}^{{\rm{LAT}}},$ respectively. ${{\boldsymbol{x}}}_{0}$ is the initial solution; the subscript zero indicates its original value. The variables were subjected to linear constraints that characterised the deformation in AP and Lat. views; equations (6) to (8) are the control vertices' movement constraints in $\vec{u},$ $\vec{v}$ and $\vec{w}$ directions, respectively. Bounds were also imposed such that the movement of the control vertices were limited to a radius of 50 mm from their initial positions (equation (9)).

Since a decent initial solution could be provided, and that the objective function was differentiable but not continuous, only local derivative-free optimisation algorithms could be implemented. Because most gradient-free local optimisation algorithms only accept bound-constraints, to resolve this issue, the augmented Lagrangian method (AUGLAG) was implemented to solve the problem, since it allowed bounds and constraints to be combined with them. Six local gradient-free optimisation algorithms were identified: the Nelder-Mead method (Nelder and Mead 1965), the PRAXIS method (Brent 1972), the Subplex method (Rowan 1990), the COBYLA algorithm (Powell 1994), the NEWUOA algorithm (Powell 2006), and the BOBYQA algorithm (Powell 2009). They were used as the sub-optimiser to evaluate the proposed deformation method, and the outcomes were compared.

2.5. Evaluation method

The proposed deformation method was evaluated using six different local gradient-free optimisation algorithms. Six sets of deformed lumbar vertebrae models, with each set consisting L1 to L5, were 3D-printed with polylactide (PLA) material, as it can be captured under x-ray. The radiographs and CT imaging of the 3D-printed models were obtained from a GE Fluorostar 7900 C-arm (General Electric, Boston, MA) and a GE Optimal CT660 (General Electric, Boston, MA), respectively.

The evaluation was carried out by comparing the mean errors of the reconstructed vertebrae using the optimisation algorithms; mean errors of the different anatomies of a vertebra were also compared. The mean error of each vertebra was defined as the averaged point-to-surface distance between the reconstructed vertebra from our proposed method and that from the CT imaging, computed using the iterative closest point (ICP) algorithm. For each lumbar segment, the mean error, its standard deviation, the root mean square (RMS) error, the 2RMS error, and the maximum error using each optimisation algorithm were measured. The averaged value of the objective function, ${\bar{d}}_{{\rm{m}}},$ and the averaged deformation duration were also recorded.

The RMS error, ${{\epsilon }}_{{\rm{RMS}}},$ represents the differences between the predicted and the observed models

$$\begin{eqnarray}&&\begin{array}{c}{{\epsilon }}_{{\rm{RMS}}}=\sqrt{\displaystyle \sum _{i=1}^{n}\displaystyle \frac{{\left({x}_{i}-{\hat{x}}_{i}\right)}^{2}}{n}}\end{array}\end{eqnarray} \tag{ 12 }$$

It is defined as the sum of the squared differences between each vertex ${x}_{i}$ to its predicted vertex ${\hat{x}}_{i}$ over the population size $n;$ the predicted ${\hat{x}}_{i}$ in our analysis is zero. The 2RMS error is defined as two times of the RMS error of the standard deviations of the point-to-surface distances of its component vertebrae models, with $m$ being the total number of vertebrae in the 6 model sets.

$$\begin{eqnarray}&&\begin{array}{c}{{\epsilon }}_{2{\rm{RMS}}}=2\cdot \sqrt{\displaystyle \sum _{i=1}^{m}\displaystyle \frac{{{\rm{SD}}}_{i}^{2}}{m}},m=30\end{array}\end{eqnarray} \tag{ 13 }$$

The deformation outcome using the six different optimisation methods are shown in figure 5, using the L2 vertebra from a model set as an example. It can be observed that errors mostly concentrate on the processes; some can be found on the vertebral body ridgeline as well.

Zoom In Zoom Out Reset image size

From the deformation outcome using the six optimisation algorithms, all but the AUGLAG + COBYLA combination produced similar mean errors ranged between 1.5 to 1.6 mm, indicating that the proposed deformation method is consistent with most local gradient-free sub-optimiser. The AUGLAG + COBYLA combination achieved the lowest mean error (1.17 mm), but since it yielded the least deformation on the template compared to the others, it was considered unsuitable, and was eliminated from our study. The AUGLAG + BOBYQA algorithm was chosen for the proposed deformation method as it took the least computation time per vertebra (4.2 seconds), and that it yielded a lower mean error compared to the others.

The averaged errors of each lumbar segment from the six model sets computed using the AUGLAG + BOBYQA were examined. As shown in table 1, the third lumbar has the lowest mean error, and it gradually increases towards the first and fifth lumbar. The same trend can be observed on the standard deviation, RMS error, and 2RMS error as well; this happens to reflect the magnitude of the absolute values of their natural pitch angles.

Table 1. Averaged errors and deviations of each lumbar segment from the six model sets computed using the AUGLAG + BOBYQA algorithm.

				Unit: mm
Lumbar segments	Mean (SD)	${{\epsilon }}_{{\rm{RMS}}}$	${{\epsilon }}_{2{\rm{RMS}}}$	Max
L1	1.66 (2.95)	2.04	2.41	8.47
L2	1.43 (2.63)	1.78	2.14	7.03
L3	1.28 (2.49)	1.63	2.04	8.74
L4	1.38 (2.63)	1.74	2.15	10.85
L5	1.76 (3.62)	2.29	2.96	13.99

Comparisons on different vertebra anatomies were performed using the AUGLAG + BOBYQA algorithm (table 2). The pedicles have the lowest mean error of 1.15 mm, whilst the spinal process has the highest mean error of 2.10 mm; the maximum error occurred at the articular processes (13.99 mm). Overall, the vertebral processes have higher error than those of the vertebral body and the pedicles: the extremities tend to concentrate at sharp ends of the deformed vertebra.

Table 2. Errors and deviations of the deformed vertebrae and their anatomies using the augmented Lagrangian method and the AUGLAG + BOBYQA algorithm.

				Unit: mm
Vertebra anatomies	Mean (SD)	${{\epsilon }}_{{\rm{RMS}}}$	${{\epsilon }}_{2{\rm{RMS}}}$	Max
Vertebra	1.50 (6.47)	1.90	2.36	13.99
Body	1.42 (5.79)	1.75	2.12	10.69
Pedicles	1.15 (4.06)	1.36	1.48	4.32
Transverse processes	1.63 (6.70)	2.02	2.45	6.94
Spinal process	2.10 (8.83)	2.62	3.23	10.85
Articular processes	1.66 (6.96)	2.05	2.54	13.99

A 3D lumbar vertebrae reconstruction method using biplanar x-rays and template models has been presented and evaluated. The AUGLAG + COBYLA combination yielded the lowest mean error, yet it did not produce any deformation to the template vertebrae; because their uniform geometries produced little regions with significant extremities, so that the error has no substantial fluctuation to increase the mean. Comparisons of the mean errors between segments of the lumbar vertebrae showed that the magnitude of the mean error is directly proportional to the absolute magnitude of the segment's natural pitch angle. During rotation estimation from the AP x-ray, the roll and yaw angles were estimated from vertebrae tilted with unknown pitch angles, thus the estimation could be inaccurate. By imposing inaccurate angle estimations to the template vertebrae and then deforming them, their pitch would not necessarily reflect what is observed on the polylactide models.

The outcomes of the proposed deformation method were compared to literatures on vertebrae reconstruction from biplanar radiographs: the template model-based NSCP method on cadaveric lumbar vertebrae (Mitulescu et al 2001) and on scoliotic thoracic and lumbar vertebrae (Mitulescu et al 2002); and statistical-based techniques by Kadoury et al (2009) and Glaser et al (2012), both were verified with scoliotic spines; the latter was the commercially available EOS system. Vertebrae models that do not contain deformities were reconstructed using the proposed method to be compared to Mitulescu et al (2001). The outcomes are shown in table 3.

For reconstructing non-deformed lumbar vertebrae, a slightly higher mean error is observed from our proposed method by about 0.17 mm, but the 2RMS error is smaller by 1 mm; this indicates that the proposed method tends to yield more local deformities.

For scoliotic cases, the mean error from our proposed method is higher than those from the statistical-based methods. Although they are capable of reconstructing accurate vertebrae models, they required extensive training beforehand, whereas the proposed method does not require any model training at all. The comparison to Mitulescu et al (2002) showed that our proposed method and the NSCP method on scoliotic cases yielded a similar mean error of 1.5 mm, and that the proposed method produced a lower 2RMS error and maximum error. The NSCP method described in Mitton et al (2000) and Mitulescu et al (2001, 2002) used 6 stereo-corresponding landmarks per vertebra to establish its orientation and a minimum of 19 non-stereo-corresponding landmarks to describe local anatomical features; the proposed method only requires 16 non-stereo-corresponding landmarks per vertebra to achieve a similar reconstruction accuracy.

The advantage of our proposed method is that radiation exposure is significantly reduced during vertebrae reconstruction compared to CT scan, as only two x-rays were taken. Although our proposed method is intended to be use intra-operatively, the transfer of files from the C-arm and the manual identification of landmarks can prolong operation time. In addition, the mean errors of the vertebral body and the pedicles in the proposed method were 1.42 mm and 1.15 mm, respectively, which are acceptable for the maximum error allowance of 2 mm for orthopaedic surgeries in the lumbar region.

Conventional orthopaedic surgeries can be performed using only C-arm images for guiding; the proposed reconstruction method would be a helpful addition to the operation when the surgeon requires a visualisation on the planned trajectories in relation to the reconstructed vertebrae models. The integration of the proposed reconstruction method to fluoroscopic-based navigation systems, such as those proposed by Lin (2019) and Chuang (2020), would be ideal for surgeons to take advantage of the convenience the intra-operative x-rays provide in tool guidance, direct visualisation, and 3D navigation altogether.

Manual identification of the anatomical landmarks on the x-rays is the main contributor to the error; poor image quality hindered accurate placement of the landmarks, hence repeatability could not be ensured. The lack of constraint to preserve the original shape of the model has led to local deformities that increased the geometrical variation. Also, the formulation of the objective function assumed orthogonal projection when measuring distances between template landmarks and their targets in AP and Lat. views, even though their depths were not the same; but it simplified the formulation, and the mean error was still comparable to literatures.

The developed programme that carried out the reconstruction could identify individual landmarks for a single vertebra based on their relative locations, and the positioning block captured in the x-rays would provide orientational information as well. The programme would inspect the selected AP- and Lat.-view landmarks for each lumbar vertebra, and only the segments with landmarks selected from both views will be used in template deformation; the user would be notified on any mis-matched landmarks before leaving the landmark-selection stage. The programme could not, however, detect whether the landmarks for a selected lumbar segment are being placed on the correct vertebra on the image; it depends on the user to identify the correct vertebra on the x-rays. In such case, the mis-matched landmarks would be used in template deformation, resulting in spurious models.

A 3D lumbar vertebrae reconstruction method using biplanar x-rays and template models has been proposed in this study. It was evaluated with six sets of 3D-printed deformed lumbar vertebrae models using the augmented Lagrangian method paired with each of the six local gradient-free optimisation algorithms. It yielded mean errors of 1.27 mm and 1.5 mm in in vitro experiments for normal and pathological lumbar vertebrae, respectively; the augmented Lagrangian method paired with the BOBYQA algorithm was deemed the most suited to the proposed deformation problem because it took the least deformation duration of 4.2 seconds, and it produced a lower error of all. Comparisons to literatures showed that our proposed method is inferior to statistical-based methods, but fairs equally with the NSCP method on clinical scoliotic cases, indicating that our proposed method is a viable alternative 3D lumbar vertebrae reconstruction method with an acceptable error, and is suitable for operations where the reconstruction error is within the allowable surgical error, and that the patient must be fully anaesthetised to prevent movements. The proposed method can be integrated to fluoroscopic-based navigation systems for further experiments in the future.

The authors would like to thank Dr. Yuan-Fu Liu and Ms Shiao-Shian Shi at NCKU Hospital and Dr Chung-Chen Shu at the Tainan Municipal An-Nan Hospital for providing access to medical equipment.

The data that support the findings of this study are available upon reasonable request from the authors.

No external funding was received for conducting this study.

The authors have no conflicts of interest to declare.