key: cord-0045038-86mbt5xi
authors: Drobny, David; Ranzini, Marta; Isaac, Amanda; Vercauteren, Tom; Ourselin, Sébastien; Choi, David; Modat, Marc
title: Towards Automated Spine Mobility Quantification: A Locally Rigid CT to X-ray Registration Framework
date: 2020-05-13
journal: Biomedical Image Registration
DOI: 10.1007/978-3-030-50120-4_7
sha: 37d54ab1851d53b6a80e447184093254e8947abb
doc_id: 45038
cord_uid: 86mbt5xi

Different pathologies of the vertebral column, such as scoliosis, require quantification of the mobility of individual vertebrae or of curves of the spine for treatment planning. Without the necessary mobility, vertebrae can not be safely re-positioned and fused. The current clinical workflow consists of radiologists or surgeons estimating angular differences of neighbouring vertebrae from different x-ray images. This procedure is time consuming and prone to inaccuracy. The proposed method automates this quantification by deforming a CT image in a physiologically reasonable way and matching it to the x-ray images of interest. We propose a proof of concept evaluation on synthetic data. The automatic and quantitative analysis enables reproducible results independent of the investigator. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this chapter (10.1007/978-3-030-50120-4_7) contains supplementary material, which is available to authorized users.

Spine mobility quantification (SMQ) describes the measurement of angles between vertebrae and their change between different positions, in order to evaluate the mobility of individual vertebrae. It is, for example, performed for scoliosis patients to determine whether vertebrae have the required mobility to be Electronic supplementary material The online version of this chapter (https:// doi.org/10.1007/978-3-030-50120-4 7) contains supplementary material, which is available to authorized users. realigned and fused to improve the curvature of the spine. SMQ is typically evaluated manually based on x-ray images. The clinician measures the angles between neighbouring vertebrae in x-ray images acquired in two opposite and extreme positions, for example flexion and extension. The change of the angular difference from one position to the other informs about the mobility of those vertebrae. Manual measurements are subjective and prone to errors, especially when a 3D movement is assessed by rotations within one plane only.

In many clinical cases a CT image is available which could be used for 3D-2D registration. This could help automating SMQ by aligning the CT image to the X-ray target pose and deriving the SMQ values from the applied deformation.

Registration of CT volumetric images to x-ray projection images is a relevant task also for several other clinical applications, such as the automated localisation of an interventional image with respect to the pre-operative CT data. Specifically, this approach can be used in image-guided spine surgery where the target vertebra has to be identified reliably [3] . According to Mody et al. more than 0.3% of spine surgery procedures have been affected by wrong-level errors which can have a severe negative impact on the patient [5] . This demonstrates the difficulty of vertebra identification even for experienced surgeons. After robust identification, a precise registration further enables navigated surgery which typically makes use of intra-operative CT (iCT) imaging. A reliable and accurate 3D-2D registration for interventional x-ray images could thus decrease errors in surgical procedures, improve accuracy of image-guided interventions or otherwise maintain accuracy while decreasing radiation exposure to patients and clinicians, compared to using iCT imaging.

Previous work in volume-projection registration (of the spine) predominantly considers rigid registration frameworks [3, 4, 8, 9, 11] . This reduces the complexity of this ill-posed problem but limits the potential accuracy and applications, as it cannot capture local deformations of soft tissue related to different positions of the subject. Non-rigid advances have been mostly focused on other body parts, for example, on lung movement [12] , and might be unrealistic in rigid anatomical structures such as bones.

To overcome the limitations of existing methods used for spine applications, we suggest a locally rigid registration framework based on the works of Arsigny et al. [1] . The proposed framework extends the poly-rigid transformation model to ensure the preservation of local rigidity during articulated movement and relies on a novel regulariser to enforce physiologically reasonable transformations. It enables better matching of images acquired with different patient position, e.g. flexion and extension, or pre-operative supine and intra-operative prone. Furthermore, the parametrisation of our model can be directly used to quantify relative positions of the vertebrae, for example for SMQ.

We demonstrate a proof of concept for spine mobility quantification based on simulated data using our registration framework. Rotation angles can be determined with high precision which motivates further development of this approach.

The proposed registration framework enables the recovery of spine movement via 3D-2D registration and the quantification of differential movement between neighbouring vertebrae. The underlying transformation model uses the polyrigid approach which combines n local sets of rigid parameters l = {l 1 , .., l n ) into a single smooth deformation field [1] . It requires a set of non-overlapping regions R i which, in our case, are given by a binary mask for each vertebra i. Based on the regions, individual weight maps w i are computed that define the local influence of each rigid transformation T i .

The registration algorithm finds an optimal choice of parameters p = (g, l), with global affine transformation parameters g and local rigid transformation parameters l, that describes the relationship between the reference x-ray image X and the CT volume V via the transformation T . The transformation T combines the poly-rigid transformation T pr (x, l) with a global affine transformation A(x, g):

To enable the comparison of 3D and 2D images, a projection P is used:

The distance measure D describes the dissimilarity of the images and is combined with a regulariser R to form the objective function F:

The parameter α controls the ratio between both terms and thus adjusts the influence of the regularisation. Minimisation of D makes the appearance of the images more similar while R promotes a anatomically realistic transformation. The numerical optimisation of F is performed iteratively via a gradient descent solver. In each step a subset of transformation parameters is evaluated, global rigid followed by local rigid parameters for each region, one at a time.

Optimising only a subset of parameters reduces the complexity of the partial derivatives and empirically leads to a more stable behaviour of the model.

For increased robustness to local minima and to increase the capture range, a multi-resolution approach is applied that performs the optimisation on multiple levels -from coarse to fine. The framework also enables 3D-3D registration of a volume V to a reference volume V ref . In this case the projection step is omitted and the objective function is therefore:

The poly-rigid transformation model is based on a set of local rigid transformation parameters l i with associated weight maps w i defining its local influence for each region i. To combine the different local rigid transformations T i into a single one, an ad hoc solution is the weighted average of the transformations' displacement fields. This does not preserve local behaviours and can lead to folding, i.e. breaking of the topology, which is to be avoided in most medical image processing. Arsigny et al. proposed a method that provides a combined transformation that is invertible and thus does not lead to folding [1] . The local rigid transformations T i are by definition homomorphisms, i.e. invertible. The average of infinitesimal small downscaled homomorphisms is also invertible as the displacement gets close to zero. The composition of homomorphisms is also a homomorphism so that the full scale transformation recovered by composition of the infinitesimal transformations is a homomorphism as well. The transformation T (t, x) is parametrised as a velocity field and integrated over time t ∈ [0, 1]. At time t = 0 the transformation is the identity and at t = 1 it is T pr .

The weighted average of the downscaled transformations is computed as

where T 1 m describes the m-th matrix square root of the rigid transformation matrix T . m has to be chosen such that the downscaled velocity field is close enough to zero [1] . This average function needs to be upscaled to get the wanted poly-rigid transformation T pr :

Arsigny et al. use an efficient scaling-and-squaring scheme to integrate the final transformation [1] . However, to maintain the desired rigid properties for each region, the weights need to be updated after each integration step as suggested by Porras et al. [6] . Therefore we use the following Euler integration scheme:

The superscript · (j) indicates that the corresponding term is updated in step j. 

Poly-rigid transformation models require weight maps w i which define the spatial influence of multiple transformations T i . Our model's weight computation is based on non-overlapping regions R i , that represent a spine segmentation.

In the region based weight computation of Porras et al. [6] , each weight w approximates the following values: w(x) = 1 if x is within the region, 0 < w(x) < 1 otherwise. This implies that each region has a global influence and the rigid properties of a region are not preserved anywhere. To guarantee the rigid transformation of a region, weights with following properties are required.

The computation of the weight map can be split into three steps. First, creating a smooth weight map that is 1 within a region and decreases outside:

where EDT denotes the Euclidean Distance Transform of the binary mask of region R i in mm (i.e. the distance to the nearest voxel of this region). The slope constant s controls how quickly the weight decays and is fixed for all regions throughout the computations.

To set the influence of all other regions R k to 0 within region R i , we multiply each region's weight by the complementary weight of all other regions:

Normalising the weights w b i to recover the maximal value of 1 within each region yields the final weight maps w i that fulfil Eq. (8):

.

(11) Figure 1 shows example weight maps computed for a CT slice containing the lowest lumbar vertebrae. This example visualises that within a region (i.e. a vertebra) the weight is one, smoothly decaying when moving away from it, and zero in each other region.

In order to guide the behaviour of the transformation model in a way that is closer to anatomical spine movement, we introduce a set of regularisation terms. These soft constraints penalise parameters that are in unexpected value ranges. The regulariser must have little influence while the parameters are within the expected value range and steeply increase when outside. Therefore the general form of each soft constraint C used is:

where b is the absolute upper bound for the value q and the exponent c controls how steep the penalty term increases. Firstly, the parameters of the global affine registration g are limited as we assume the general direction of x-ray acquisition is known and thus, the volume should not translate or rotate too much away from the expected view.

Secondly, the local rigid parameters l are constrained via the parameter difference of pairs of neighbouring vertebrae, separate for rotation and translation:

This way neighbouring vertebrae cannot translate or rotate too far from their initial relative position. This constraint guarantees that the spine moves consistently and retains its integrity. An unregularised optimisation step might lead to local rigid parameters that would move regions into another. The overall poly-rigid transformation avoids this kind of folding but compromises the rigidity of the underlying transformations to achieve this. The third regulariser R o is used to discourage such optimisation steps by penalising transformations that cause voxel overlap in order to achieve diffeomorphic transformations while preserving the local rigidity.

The sum of all individual terms leads to the overall regulariser R as in Eq. (3):

Data from the Spineweb database which consist of CT images and corresponding vertebra segmentation have been used for the experiments described in this work [13] . A CT image is cropped to the lumbar area including six vertebrae above the sacrum and resampled to a 1 mm isotropic image of 162 × 162 × 312 voxel. The rigid parameters of the vertebrae are user-defined to simulate two movements: flexion and extension. The original CT image is considered as the starting position and angles are interpreted as differential angles to this position. The image axes x, y, and z correspond to the lateral, anterior-posterior, and superior-inferior patient axis, respectively. X-ray projections are simulated from the deformed CT images in lateral direction (along the x-axis) and are used as the reference images X while the initial CT volume is then used as moving image V. The projection method P used for the synthetic experiments is a basic averaging of intensity values along a given direction, which can be easily replaced by other approaches to account for more clinically realistic scenarios, like those presented by Unberath et al. [7] . Table 1 shows the rotation angles chosen for each vertebra in the flexion and extension simulation. For this experiment only rotations around the x-and yaxis are considered, as those are the angles clinically measured for SMQ and rotations around z-axis are minor in the lumbar spine. The values are in the typical range observed in clinical practice, for example as reported by Wilke et al. [10] . Figure 2 shows the flexion, extension and initial position of the synthetic data as 3D rendered visualisations. The vertebrae are deformed rigidly while the surrounding soft tissue, like the arteries, are non-linearly deformable.

In our experiment, images in two different poses are target of the registration: flexion and extension. This resembles the application of SMQ where usually two opposite extremal positions are compared. For each pose, two experiments are evaluated: (1) 3D-3D registration-matching the original CT V to the synthetic CT image V ref and (2) 3D-2D registration-matching the CT image V to the synthetic x-ray image X . Two metrics are used to evaluate the registration results: the error of the estimated angle difference and the average error of the recovered displacement field at the location of the vertebrae. The displacement field error gives information about whether the vertebrae are in the right position while the angle errors describe how well the orientation of the vertebrae was recovered, which is the focus of SMQ.

For performance reasons, the poly-rigid transformation computed by Eq. (7) is used only as a post processing step and we use the direct weighted average during the optimisation, as was suggested by Commowick et al. [2] . We used the mean of squared difference as similarity measure D during our experiments. The parameters are either motivated by clinical movement ranges or chosen empirically: four levels of multi-resolution approach, with the finest resolution of 1 mm, the weight constant s = 0.05, the constraint scaling parameter α = 100. Referring to Eqs. 12 and 16, the first two constraint functions have the exponent c = 6, scaling γ = 1 and the boundaries b of 20 mm, 30 • for the global rigid 

Our experiments show that with both 3D-3D as well as 3D-2D registration, the parameters of the synthetic images can be recovered. In Table 2 the quantitative results of all experiments are summarised. Using the 3D-3D registration higher accuracy was achieved compared to the 3D-2D cases which is due to the higher information retained in a volume compared to its projection. In the 3D-3D experiments, the average of the mean absolute angle error was reduced from the initial 2.47 • to 0.12 • and from 2.70 • to 0.16 • for extension and flexion respectively. In the 3D-2D case, the respective average errors are 0.59 • and 0.42 • . The deformation field error (DFE) is also reported as a measure of registration accuracy in Table 2 . The DFE confirms the higher accuracy of the 3D-3D approach compared to the 3D-2D registration.

For the spine mobility quantification the angular difference between the two positions for each vertebra is of interest. This can be directly derived from the individual registration results as the difference of the respective angles. 

The non-rigid alignment of a volumetric to a projection image is challenging because several degrees of freedom affecting the 3D volume have to be recovered by only comparing 2D images. We presented a framework based on a locally rigid transformation model that enables 3D-3D and 3D-2D registration of spine images. These synthetic experiments show that both 3D-3D and 3D-2D registration can be recovered, demonstrating that the registration with our transformation model can lead to clinically useful results for spine mobility quantification. This framework offers the full processing pipeline needed for different applications and is flexible to adaptations required by the use of real clinical data. We also introduced simple constraints to guide a consistent movement of the spine while optimising individual vertebrae parameters iteratively. Our experiments are limited by only considering synthetic x-ray images generated from deformed CT scans. Clinical x-ray images are more challenging to register and thus need specialised image similarity measures. Future work will thus focus on the identification of a measure suitable for clinical data, and on the benchmarking of the proposed approach against the current clinical practice on real patient data. Furthermore we will use a higher level parametrisation of the spine shape, e.g. using principal components of a spine shape model, to optimise the position of multiple vertebrae simultaneously in a consistent way. This will be helpful to capture larger deformations especially in an extended field of view. Such a model could also be used to extract population statistics like the principal modes of spine variability and movement patterns. As the proposed framework proved effective also for 3D-3D registration, it will be tested on further clinical applications, such as spine CT-MRI registration. Finally, as this pipeline requires vertebra segmentations, automated segmentation techniques (e.g. using deep learning approaches) will be explored to complete the registration framework.

A fast and log-euclidean polyaffine framework for locally linear registration

An efficient locally affine framework for the smooth registration of anatomical structures

3D-2D image registration for target localization in spine surgery: investigation of similarity metrics providing robustness to content mismatch

A versatile intensity-based 3D/2D rigid registration compatible with mobile C-arm for endovascular treatment of abdominal aortic aneurysm

The prevalence of wrong level surgery among spine surgeons

Locally affine diffeomorphic surface registration and its application to surgical planning of fronto-orbital advancement

DeepDRR -a catalyst for machine learning in fluoroscopyguided procedures

Fully automated 2D-3D registration and verification

A comparison of two novel similarity measures based on mutual information in 2D/3D image registration

The role of prosthesis design on segmental biomechanics semi-constrained versus unconstrained prostheses and anterior versus posterior centre of rotation

Evaluation of similarity measures for use in the intensity-based rigid 2D-3D registration for patient positioning in radiotherapy

Thorax x-ray and CT interventional dataset for nonrigid 2D/3D image registration evaluation

Detection of vertebral body fractures based on cortical shell unwrapping