Keywords

1 Introduction

The main objective in orthopedic reduction surgery is to restore the correct alignment of the dislocated or fractured bone. In both unilateral and bilateral fractures, surgeons attempt to re-align the fractures to their natural biological alignment. In the majority of cases, there are no available anatomical imaging data prior to injury, and CT scans are only acquired after the patient is injured to identify the fracture type and plan the intervention. Therefore, no reference exists to identify the correct and natural alignment of the bone fragments. Instead, surgeons use the opposite healthy side of the patient as reference, and aim at producing symmetry across the sagittal plane [1]. It is important to mention that, although the healthy pelvic bone is not entirely symmetric, surgeons aim at aligning the bone fragments to achieve structural and functional symmetry. For orthopedic traumatologists, the correct length, alignment, and rotation of extremities are also verified by comparing to the contralateral side. Examples of other fields of surgery that use symmetry for guidance include crainiofacial [2] and breast reconstruction surgeries [3].

Self-symmetry assessment is only achievable if the fractures are not bilateral, so that the contralateral side of the pelvis is intact. According to pelvis fracture classification, a large number of fracture reduction cases are only unilateral [4]. Consequently, direct comparison of bony structures across the sagittal plane is possible in a large number of orthopedic trauma interventions.

CT-based statistical models from a population of data are particularly important when patient-specific pre-operative CT images are not present. In this situation, statistical modeling and deformable registration enable 3D understanding of the underlying anatomy using only 2D intra-operative imaging. Statistical shape models are used to extrapolate and predict the unknown anatomy in partial and incomplete medical images [5]. In the aforementioned methods, instead of patient self-correlation, relations to a population of data are exploited for identifying missing anatomical details.

In this work, we hypothesize that there is a high structural correlation across the sagittal plane of the pelvis. Quantitative 3D measurements on healthy pelvis data indicate that \(78.9 \%\) of the distinguishable anatomical landmarks on the pelvis are symmetric [6], and the asymmetry in the remaining landmarks are still tolerated for fracture reduction surgery [7]. To exploit the partial symmetry, we automatically detect the desired plane of symmetry using Tukey’s biweight distance measure. In addition, a novel regularization term is designed that ensures a similar distribution of bone density on both sides of the plane. Regularization is important when the amount of bone dislocation is large, and Tukey’s cost cannot solely drive the symmetry plane to the optimal pose. After identifying the partial symmetry, the CT volume is mirrored across the symmetry plane which then allows simulating the ideal bone fragment configurations. This information is provided intra-operatively to the surgeon, by overlaying the C-arm X-ray image with a forward-projection of the mirrored volume.

The proposed approach relies on pre-operative CT scans of the patient. It is important to note that acquiring pre-operative CT scans is standard practice in severe trauma and fracture reduction cases. Therefore, it is valid to assume that pre-operative imaging is available for the types of fractures discussed in this manuscript, namely illiac wing, superior and inferior pubic ramus, lateral compression, and vertical shear fractures. In this paper, we introduce an approach that enables the surgeon to use patient CT scans intra-operatively, without explicitly visualizing the 3D data, but instead using 2D image augmentation on commonly used X-ray images.

2 Materials and Methods

2.1 Problem Formulation

The plane of symmetry of an object \(O \subseteq \mathbb {R}^3\) is represented using an involutive isometric transformation \(\mathbf {M}(g) \in E_o(3)\), such that \(E_o(3) = \{ h \in E(3) \vert h (O) = O \}\), where E(3) is the 3D Euclidean group, consisting of all isometries of \(\mathbb {R}^3\) which map \(\mathbb {R}^3\) onto \(\mathbb {R}^3\). The transformation \(\mathbf {M}(g)\) mirrors the object O across the plane as; \(\mathbf {o_{-x}} = \mathbf {M}(g) \, \mathbf {o_{x}}\), where \(\mathbf {o_x} \subseteq \mathbb {P}^3\) and \(\mathbf {o_{-x}} \subseteq \mathbb {P}^3\) are sub-volumes of object O on opposite sides of the symmetry plane, and are defined in the 3D projective space \(\mathbb {P}^3\). Assuming the plane of symmetry is the Y-Z plane, \(\mathbf {M}(g)\) is given by \(\mathbf {M}(g) := g \, \mathbf {F}_x \, g^{-1}\), where \(\mathbf {F}_x\) is reflection about the X-axis, and \(g \in SE(3)\), where SE(3) is the 3D special Euclidean group. We propose to estimate \(\mathbf {M}(g)\) by minimizing a distance function \(D (\mathbf {M}(g))\) as:

$$\begin{aligned} \mathop {{{\mathrm{arg\,min}}}}\limits _{g} D(\mathbf {M}(g)) := d_{I}(\mathbf {o_{x}}, \mathbf {M}(g)\mathbf {o_{x}}) + \lambda \; d_{D}(\mathbf {o_{x}}, \mathbf {M}(g)\mathbf {o_{x}}), \end{aligned}$$
(1)

combining an intensity-based distance measure \(d_{I}(.)\), and a regularization term \(d_{D}(.)\) based on the bone density distribution. The term \(\lambda \) in Eq. 1 is a relaxation factor, and \(\lambda , d_{I}(.), d_{D}(.) \in \mathbb {R}\). Derivations of \(d_{I}(.)\) and \(d_{D}(.)\) are explained in Sects. 2.2 and 2.3, respectively. In Sect. 2.4, we suggest an approach to incorporate the knowledge from the plane of symmetry, and provide patient-specific image augmentation in fracture reduction interventions.

2.2 Robust Loss for Estimation of Partial Symmetry

The CT data of a pelvis only exhibits partial symmetry, as several regions within the CT volume may not have a symmetric counterpart on the contralateral side. These outlier regions occur either due to dislocation of the bone fragments, or asymmetry in the natural anatomy. To estimate the plane of “partial symmetry”, we suggest to minimize a disparity function that is robust with respect to outliers, and only considers the partial symmetry present in the volumetric data. The robustness to outlier is achieved by down-weighting the error measurements associated to potential outlier regions. To this end, we estimate the plane of partial symmetry by minimizing Tukey’s beweight loss function defined as [8]:

$$\begin{aligned} \begin{aligned}&d_{I}(\mathbf {o_{x}}, \mathbf {M}(g) \mathbf {o_{x}}) = \sum _{i=1}^{ |\Omega _s|} \frac{\rho (e_i(\mathbf {M}(g)))}{|\Omega _s|},\\&\rho (e_i(\mathbf {M}(g))) = {\left\{ \begin{array}{ll} e_i(\mathbf {M}(g))\left[ 1 - \left( \frac{e_i(\mathbf {M}(g))}{c} \right) ^2 \right] ^2 &{}; |e_i(\mathbf {M}(g))| \leqslant c,\\ 0 &{}; \text {otherwise}, \end{array}\right. } \end{aligned} \end{aligned}$$
(2)

with \(\Omega _s\) being the spatial domain of CT elements. The threshold of assigning data elements as outlier is defined by a constant factor c that is inversely proportional to the down-weight assigned to outliers. As suggested in the literature, \(c = 4.685\) provides high asymptotic efficiency [8]. The residual error for the i-th voxel element is \(e_i(\mathbf {M}(g))\) and is defined as following:

$$\begin{aligned} e_i(\mathbf {M}(g)) = \frac{| \text {CT}(o_{x_i}) - \text {CT}(\mathbf {M}(g) o_{x_i}) |}{S}. \end{aligned}$$
(3)

In Eq. 3, \(o_{x_i}\) is the i-th voxel element, \(\text {CT}(o_{x_i})\) is the intensity of \(o_{x_i}\), and S is the scaling factor corresponding to the median absolute deviation of residuals.

2.3 Bone Density Histogram Regularization

In fractured bones, the distribution of bone material inside the body will remain nearly unaffected. Based on this fact, we hypothesize that the distribution of bone intensities, i.e. histograms of bone Hounsfield Unit (HU), on the opposite sides of the plane of symmetry remains similar in presence of fracture (example shown in Fig. 1). Therefore, we design a regularization term based on normalized mutual information as follows [9]:

$$\begin{aligned} d_{D}(\mathbf {o_{x}}, \mathbf {M}(g) \mathbf {o_{x}}) = - \frac{H(\text {CT}(\mathbf {o_{x}})) + H(\text {CT}(\mathbf {M}(g) \mathbf {o_{x}}))}{H(\text {CT}(\mathbf {o_{x}}), \text {CT}(\mathbf {M}(g) \mathbf {o_{x}})))}, \end{aligned}$$
(4)

where H(.) is the entropy of voxels’ HU distribution. Minimizing the distance function in Eq. 4 is equivalent to increasing the similarity between the distributions of bone on the opposing sides of the plane of partial symmetry.

Fig. 1.
figure 1

Bone intensity histograms are shown in (a–b) and (c–d) before and after estimating the plane of partial symmetry.

2.4 Patient-Specific Image Augmentation

After estimating the plane of partial symmetry, the CT volume is mirrored across this plane to construct a patient-specific CT volume representing the bony structures “as if they were repaired”. It is important to note that, although the human pelvic skeleton is not entirely symmetric, it is common in trauma interventions to consider it as symmetric, and use the contralateral side as reference.

Fig. 2.
figure 2

Workflow for patient-specific image augmentation based on partial symmetry

To assist the orthopedic surgeon in re-aligning bone fragments, we propose to augment intra-operative X-ray images with the bone contours from the mirrored CT volume. This step requires generation of digitally reconstructed radiographs (DRRs) from views identical to the one acquired intra-operatively using a C-arm. Hence, 2D/3D image registration is employed to estimate the projection geometry between the X-ray image and CT data. This projective transformation is then used to forward-project and generate DRR images from the same viewing angle that the X-ray image was acquired. Finally, we augment the X-ray image with the edge-map acquired from the DRR that will then serve as a road-map for re-aligning the fragmented bones. The proposed workflow is shown in Fig. 2.

3 Experimental Validation and Results

We conducted experiments on synthetic and real CT images of healthy and fractured data. For all experiments, the optimization was performed using bound constrained by quadratic approximation algorithm, where the maximum number of iterations was set to 100. The regularization term \(\lambda \) in Eq. 1 was set to 0.5 which allowed \(d_I(.)\) to be the dominant term driving the similarity cost, and \(d_{D}(.)\) to serve as a data fidelity term.

Performance Evaluation on Synthetic Data: A synthetic 3D data of size \(100^3\) voxels and known plane of symmetry was generated. The synthetic volume was perfectly symmetric along the Y-Z plane, where the intensities of each voxel was proportional to its Y and Z coordinate in the volume. We evaluated the performance of the Tukey-based cost \(d_I(.)\) with respect to noise and outliers, and compared the outcome to NCC-based cost. The results are shown in Fig. 3.

Fig. 3.
figure 3

The performance of Tukey-based distance measure is evaluated against noise level and amount of outlier. The horizontal axis indicates the relative pose of the symmetry plane at the initial step with respect to the ground-truth, where the initialization parameters were increased with the increments of 2 voxels translation and \(2^{\circ }\) rotation along each 3D axis. The amount of outlier was varied between \(0\%\) to \(30\%\) of the size of the entire volume, and the Gaussian noise between \(0 \%\) to \(25\%\) of the highest intensity in the volume. (a–d) and (e–h) show the performance of NCC and Tukey robust estimator, respectively. It is important to note that different heat-map color scales are used for NCC and Tukey to demonstrate the changes within each sub-plot.

Estimating the Plane of Partial Symmetry on non-Fractured Data: The plane of partial symmetry was estimated for twelve CT datasets with no signs of fractures or damaged bone. Four of the volumes were lower torso cadaver CT data, and eight were from subjects with Sarcoma. After estimating the plane of partial symmetry, the CT volumes were mirrored across the estimated plane. We then identified the following 4 anatomical landmarks and measured the distance between each landmark on the original volume to the corresponding landmark on the mirrored CT volume: L\(_1\): anterior superior iliac spine, L\(_2\): posterior superior iliac spine, L\(_3\): ischial spine, and L\(_4\): ischial ramus. Results of this experiment using NCC, Tukey robust estimator \(d_I(.)\), and regularized Tukey \(d_I(.) + \lambda d_D(.)\) distance functions are shown in Table 1.

Table 1. Errors in estimation of partial symmetry were measured using four anatomical landmarks. The values in the table are in mm units, and are shown as mean ± SD. The final row are the results of our proposed method.

Estimating the Plane of Partial Symmetry on Fractured Data: We simulated three different fractures - i.e. iliac wing, pelvic ring, and vertical shear fractures (shown in Fig. 4(a–c)) - and evaluated the performance of the proposed solution in presence of bone dislocation. These three fractures were applied to three different volumes, and in total nine fractured CT volumes were generated. The error measurements are presented in Table 2.

Table 2. Error measurements on fractured data in mm units.
Fig. 4.
figure 4

(a–c) are simulated iliac wing, pelvic ring, and vertical shear fractures, respectively. Image augmentations are shown in (d–f). The red arrows indicate the fracture location in the pelvis. The green contours represent the desired bone contours to achieve bilateral symmetry. (g–h) show the 3D visualization of the fractures, where red color indicates regions that were rejected as outliers using Tukey’s robust estimator. The intra-operative X-ray image in (j) is augmented with the edge-map of the DRR generated from the mirrored volume using the projection geometry estimated from 2D/3D X-ray to CT image registration(k).

Intra-operative Image Augmentation: After estimating the plane of partial symmetry, we mirrored the healthy side of the pelvis across the plane of partial symmetry. In Fig. 4(d–f) we present the superimposition of the fractured data shown in Fig. 4(a–c) with the edge-map extracted from gradient-weighted DRRs. Moreover, we preformed 2D/3D intensity-based registration between the pre-operative CT and the intra-operative X-ray (Fig. 4j), and used the estimated projection geometry to generate DRRs and augment the X-ray image of the fractured bone. The augmentation is shown in Fig. 4k.

4 Discussion and Conclusion

This work presents a novel method to estimate partial structural symmetry in CT images of fractured pelvises. We used Tukey’s biweight robust estimator which prevents outlier voxel elements from having large effects on the similarity measure. Moreover, Tukey’s distance function is regularized by enforcing high similarity in bilateral bone HU distribution. The experimental results on synthetic data indicate that Tukey-based similarity cost outperformed NCC-based similarity cost substantially in the presence of noise and outliers. The results in Table 1 show an average landmark error of 5.78 mm and 3.46 mm using Tukey- and regularized Tukey-based cost, respectively. Similarly for fractured data presented in Table 2, the mean error reduced from 7.91 mm to 3.50 mm after including the regularization term.

In conclusion, we proposed to incorporate the knowledge from partial symmetry and provide intra-operative image augmentation to assist orthopedic surgeons in re-aligning the bone fragments with respect to bilateral symmetry. Our work relies on pre-operative CT images and is only applicable to surgical interventions where pre-operative 3D imaging exist. This solution enables patient-specific image augmentation which is not possible using statistical atlases. Using atlases for this application requires a large population of patient pelvis data for different age, sex, race, disease, etc. which are not available.