A scale-space curvature matching algorithm for the reconstruction of complex proximal humeral fractures

Medical Image Analysis 43 (2018) 142–156

Contents lists available at ScienceDirect

Medical Image Analysis

journal homepage: www.elsevier.com/locate/media

A scalespace curvature matching algorithm for the reconstruction of complex proximal humeral fractures

Lazaros Vlachopoulos a , b , ∗, Gábor Székely b , Christian Gerber c , Philipp Fürnstahl a

a Computer Assisted Research and Development Group, Balgrist University Hospital, University of Zurich, Forchstrasse 340, CH8008 Zurich, Switzerland b Computer Vision Laboratory, ETH Zurich, Sternwartstrasse 7, CH8092 Zürich, Switzerland c Department of Orthopaedics, Balgrist University Hospital, University of Zurich, Forchstrasse 340, CH8008 Zurich, Switzerland

a r t i c l e i n f o

Article history: Received 18 January 2017 Revised 26 October 2017 Accepted 26 October 2017 Available online 27 October 2017

Keywords: Scalespace Curvature Fracture reconstruction Fracture line Proximal humerus

a b s t r a c t

The optimal surgical treatment of complex fractures of the proximal humerus is controversial. It is proven that best results are obtained if an anatomical reduction of the fragments is achieved and, therefore, computerassisted methods have been proposed for the reconstruction of the fractures. However, com plex fractures of the proximal humerus are commonly accompanied with a relevant displacement of the fragments and, therefore, algorithms relying on the initial position of the fragments might fail. The state oftheart algorithm for complex fractures of the proximal humerus requires the acquisition of a CT scan of the (healthy) contralateral anatomy as a reconstruction template to address the displacement of the fragments. Poseinvariant fracture line based reconstruction algorithms have been applied successful for reassembling broken vessels in archaeology. Nevertheless, the extraction of the fracture lines and the necessary computation of their curvature are susceptible to noise and make the application of previous approaches difficult or even impossible for bone fractures close to the joints, where the cortical layer is thin. We present a novel scalespace representation of the curvature, permitting to calculate the correct alignment between bone fragments solely based on corresponding regions of the fracture lines. The frac tures of the proximal humerus are automatically reconstructed based on iterative pairwise reduction of the fragments. The validation of the presented method was performed on twelve clinical cases, surgically treated after complex proximal humeral fracture, and by cadaver experiments. The accuracy of our ap proach was compared to the stateoftheart algorithm for complex fractures of the proximal humerus. All reconstructions of the clinical cases resulted in an accurate approximation of the pretraumatic anatomy. The accuracy of the reconstructed cadaver cases outperformed the current stateoftheart algorithm. © 2017 Elsevier B.V. All rights reserved.

1. Introduction

The treatment of comminuted fractures of the proximal humerus is challenging and the optimal procedure remains con troversial ( Cvetanovich et al., 2016; Gerber et al., 2004 ). Open re duction and internal fixation using conventional or locking plates is the mainstay of therapy for the young and active patient ( Gerber et al., 2004; Grubhofer et al., 2016 ), while best results are ob tained if anatomical or near anatomical reduction can be achieved ( Gerber et al., 2004 ). Anatomical reduction is a prerequisite for a jointpreserving surgical treatment of a fractured proximal humerus. If anatomical reduction cannot be obtained, joint re

∗Corresponding author at: Computer Assisted Research and Development Group, Balgrist University Hospital, University of Zurich, Forchstrasse 340, CH8008 Zurich, Switzerland. Email addresses: lazaros.vlachopoulos@balgrist.ch (L. Vlachopoulos),

szekely@vision.ee.ethz.ch (G. Székely), christian.gerber@balgrist.ch (C. Gerber), philipp.fuernstahl@balgrist.ch (P. Fürnstahl).

placement has to be considered ( Gerber et al., 1998 ). The op tions for replacement surgery of the shoulder joint include hemi arthroplasty, anatomic total shoulder arthroplasty and reverse total shoulder arthroplasty (RTSA) ( Cuff and Pupello, 2013; Cvetanovich et al., 2016; Fucentese et al., 2014; Grubhofer et al., 2016 ) with a current trend from hemiarthroplasty towards RTSA for com plex humeral fractures in the elderly ( Cvetanovich et al., 2016; Grubhofer et al., 2016 ). The main reason of this current trend is that, despite promising initial reports of the hemiarthroplasty ( Neer, 1970 ), less satisfactory or even disappointing results have been reported ( Shukla et al., 2016 ). Current literature suggests that RTSA might result in better clinical outcomes than hemi arthroplasty, due to the decreased reliance on tuberosity healing of the RTSA ( Shukla et al., 2016 ). Nevertheless, the most important consensus across all surgical treatment options is, that the func tional outcome is better with anatomical fixation of the tuberosi ties ( Anakwenze et al., 2014; Boileau et al., 2002; Fucentese et al., 2014; Gallinet et al., 2009; Gerber et al., 2004; Grubhofer et al., 2016; Huffman et al., 2008 ). Therefore, it seems clearly justified

https://doi.org/10.1016/j.media.2017.10.006 13618415/© 2017 Elsevier B.V. All rights reserved.

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 143

that major effort should be made to achieve an anatomical reduc tion of the tuberosities. The benefits of computerassisted preoperative simulation and intraoperative navigation is well accepted in joint replacement surgery ( Iannotti et al., 2014; Levy et al., 2014; Nguyen et al., 2009 ) and for corrective osteotomies after malunited fractures of the humerus ( Murase et al., 2008; Vlachopoulos et al., 2016b ). Computerassisted approaches are promising, especially, since it is difficult or even impossible to preoperatively plan the or thopaedic procedure using only radiographs or computed tomog raphy (CT) analysis ( Vlachopoulos et al., 2016a ). In presence of a complex fracture of the proximal humerus, the ultimate goal of these approaches should be the restoration of the normal humeral anatomy. However, the fundamental prerequisite to ap ply computerassisted navigation in the surgery to fractures of the proximal humerus is the preoperative reconstruction of the frac tures ( Cui et al., 2007; Fürnstahl et al., 2012; McBride and Kimia, 20 03; Papaioannou and Theoharis, 20 03; Üçoluk and Toroslu, 1999 ). Hitherto, one method ( Fürnstahl et al., 2012 ) has been published and validated for the computerassisted reconstruction of complex proximal humerus fractures ( JimenezDelgado et al., 2016 ). Fürnstahl et al. (2012) demonstrated that their algorithm al lows accurate reconstruction of the pretraumatic anatomy. How ever, the main drawback of the method of Fürnstahl et al. (2012) is the dependency on the healthy contralateral bone model as a reconstruction template. A further computerassisted method for the treatment of complex proximal humeral fractures via hemi arthroplasty was developed and validated by cadaver experiments ( Bicknell et al., 2007 ). The alignment of the shoulder prosthe ses and the tuberosity fragment was assessed by manual mea surements of characteristic landmarks. Here, also the contralat eral anatomy was proposed as a reconstruction template for the use in a clinical setting. However, existing bilateral differences in the humeral anatomy ( DeLude et al., 2007; Vlachopoulos et al., 2016a ) or the presence of a pathological altered contralateral anatomy (e.g., after a proximal humeral fracture or a joint re placement surgery) might limit the clinical application of both methods. The task to be performed is similar to the assembly of a jigsaw puzzle as illustrated in Fig. 1 , also introducing the ter minology used throughout this paper. Poseinvariant reconstruc tion algorithms have been successfully developed in archaeology for the reassembling task of broken vessels and relicts ( McBride and Kimia, 2003; Papaioannou and Theoharis, 2003; Üçoluk and Toroslu, 1999 ). The geometric reconstruction was performed by matching individual fracture surfaces using 3D curvature matching methods. The presented method builds on this idea from the ap proaches in archaeology. However, as the data acquisition scans in clinical practice is based on CT, the noise is greater that in archae ology, where laser scanning is used. Clinical data are characterized by a limited resolution (inplane and axial resolution of 0.4 mm or worse) ( Lecouvet et al., 2008 ), resulting in partial volume ef fects and diffuse fracture lines, in contrast to laserscanned data with an isotropic resolution of 0.05 mm or less. Furthermore, bone fracture surfaces tend to be highly irregular ( Thomas et al., 2011 ). Therefore, the adoption of the archaeological approaches ( McBride and Kimia, 2003; Papaioannou and Theoharis, 2003; Üçoluk and Toroslu, 1999 ) for clinical application, i.e., bone fracture reconstruc tion, is not straightforward ( Thomas et al., 2011 ). In this paper, we present a novel method for the fully auto mated reconstruction of proximal humeral fractures, requiring only the information of the fracture surfaces. The idea is to use a scale space representation of the curvature of the corresponding frac ture lines, which permits determining the correct alignment be tween fragments. The key novelties about the proposed method are:

Fig. 1. Terminology for the present paper illustrated on a jigsaw puzzle. The frac ture surfaces (yellow surface) represent the break through the cortex of the prox imal humeral fragments and are simplified by fracture lines (red lines). The frag ment surface (blue surface) corresponds to the unfractured outer cortical layer of the fragment. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

• Bspline based fracture line representation with a tailored weighting schema for proximal humerus fractures. In con trast to previous work that investigated the reassembling of twodimensional planar fragments ( McBride and Kimia, 2003 ) or curvature matching for binary twodimensional images ( Cui et al., 2007 ) our method was investigated on 3D fracture line representation of bone fractures.

• Curvature ScaleSpace: Our scalespace representation incorpo rates the local shape of a curvature (i.e., concave or convex), in contrast to previous work ( Cui et al., 2007; McBride and Kimia, 20 03; Papaioannou and Theoharis, 20 03; Üçoluk and Toroslu, 1999 ), permitting to reduce the number of incorrect matches in scalespace and to increase the robustness of the algorithm.

• Curvature Matching Algorithm: The developed similarity crite ria detect correspondences based on the shape of curvature in scale space and the introduced normalized measure of reduc tion replaces the calculation of the torsion as a signature of a curvature – as proposed by Papaioannou and Theoharis (2003 ) and Üçoluk and Toroslu (1999 ), being robust against difference in magnitude (diffuse fracture lines, noise).

• A graphbased algorithm for robust merging of reduced bone fragments allowing automatic iterative fracture reconstruction.

• A reduction algorithm that determines the best solution based on all performed reconstructions and a warning mechanism, i.e., if only a partial fracture reconstruction is performed.

The proposed method was evaluated clinically on a consecu tive series of patients treated with proximal humerus fractures and on four artificially created fractures on cadaveric humeri. Best to our knowledge, it is the largest published set of computer reconstructed fractures of the proximal humerus. In addition, we compared our reconstruction results with the current stateofthe art algorithm ( Fürnstahl et al., 2012; JimenezDelgado et al., 2016 ). In the following, we will give a brief overview of computer assisted techniques for the simulation of fracture reduction. In Section 2 an overview is presented and the details of our approach are described. The clinical evaluation and the results of cadaver ex periments are presented in Section 3 . We discuss the method in

144 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Section 4 . Finally, Section 5 summarizes the major conclusions of the work.

1.1. Related work

An accurate preoperative assessment of fragment displacement is crucial for a successful restoration of a fracture ( Fürnstahl et al., 2012 ). However, the literature regarding computerassisted recon struction of bone fractures is relatively sparse, compared to the very large body of research on the topics of bone segmentation and medical image registration ( Thomas et al., 2011 ). Jiménez Delgado et al. (2016) recently published a comprehensive review article, summarizing current approaches for bone fracture reduc tion planning ( Albrecht and Vetter, 2012; Chowdhury et al., 2009; Fürnstahl et al., 2012; Moghari and Abolmaesumi, 2008; Okada et al., 2009; Thomas et al., 2011; Willis et al., 2007; Winkelbach and Wahl, 2008; Zhou et al., 2009 ). Some of these approaches rely on a reconstruction template, i.e. the contralateral bone model ( Bicknell et al., 2007; Fürnstahl et al., 2012; Okada et al., 2009 ) or a statistical shape model ( Albrecht and Vetter, 2012; Moghari and Abolmaesumi, 2008 ) is required to calculate the reduction. DeLude et al. (2007) and Vlachopoulos et al. (2016a) verified that the contralateral humeral anatomy might be a reliable template for some geometric charac teristics (i.e., the humeral head size and the humeral length). How ever, due to the presence of intraindividual differences, in partic ular, differences in axial torsion, Vlachopoulos et al. (2016a) con cluded that preoperative planning approaches, targeting the recon struction of complex proximal humerus fractures should not rely blindly on the contralateral anatomy. Other approaches align the fragments based on the characteris tics of the fracture surfaces ( Chowdhury et al., 2009; Kronman and Joskowicz, 2013; Okada et al., 2009; Willis et al., 2007; Winkelbach et al., 2004; Zhou et al., 2009 ). Most of the presented approaches have in common that an Iterative Closest Point (ICP) based algo rithm is used to perform the reduction of the fragments ( Jimenez Delgado et al., 2016 ). The tendency of the ICP to fall into local minima might be particularly problematic in case of the proximal humerus due to the almost spherical shape ( Fürnstahl et al., 2012 ). Further methods have been developed, that combine the use of a reconstruction template and the fractured surfaces to tackle this problem ( Fürnstahl et al., 2012; Okada et al., 2009 ). The contralat eral anatomy is used, thereby, to generate a set of initial transfor mation close to the true parameter values. Furthermore, the frag ments at the proximal humerus are often considerably displaced and malrotated and, consequently, approaches relying on the initial position of the fragments ( Kronman and Joskowicz, 2013; Okada et al., 2009 ) can likely fail. Therefore, Fürnstahl et al. (2012) , pro posed to perform first a global but coarse preregistration step that is independent of the initial pose of the fragments. Thereafter, a re finement step was applied using stateoftheart local optimization techniques. However, their global preregistration again required the use of the contralateral anatomy as a shape prior. A fracture reconstruction algorithm that does not rely on the initial pose of the fragments would be an alternative solution to tackle the problematic of the considerably displaced fragments of a proximal humerus fracture. Similar methods were developed in archaeology for the reconstruction of broken vessels ( Papaioannou and Theoharis, 2003; Üçoluk and Toroslu, 1999 ). These methods are based on the signature of a 3D curve, i.e., arc length, cur vature and torsion, as introduced by Kishon and Wolfson (1987) , and are closely related to our approach. Papaioannou and Theo haris (2003) calculated first pairwise reductions based on the fracture facets of the fragments and applied a facet boundary curve matching to reduce the search space. The best assembly was obtained by determining the set of fragment combination

that resulted in the smallest accumulative matching error. How ever, the main disadvantage of the method of Papaioannou and Theoharis (2003) is the requirement of nearly planar fracture facets. Furthermore, the limited resolution and noise of the clin ical data make the application of the archaeological approaches ( Papaioannou and Theoharis, 2003; Üçoluk and Toroslu, 1999 ) based on 3D curve matching probably more difficult ( Thomas et al., 2011 ).

2. Method

2.1. Overview of the algorithm

The overall workflow of our method consists of three modules as illustrated in Fig. 2 . In the first module, we perform the segmen tation task. The input of the first module is a CT scan of a proximal humeral fracture. Note, that the description of the segmentation task (with the annotation of the fracture surfaces) is used for the overview of the planning workflow and is not part of the devel oped method. The input of the second module are the triangular surface mod els of the cortical layer of the fragments (2.2). The characteristic of the fragments used for the fracture reconstruction are analysed as described in 3.2. In a first step, the fractured surfaces of the frag ments are converted to a different representation based on con nected line segments, herein called fracture lines as described in 2.3 . These fracture lines are the base for all subsequent steps of the fracture reconstruction algorithm. Thereafter, the scalespace rep resentation of the curvature is calculated in a local neighbourhood for each point of a fracture line in Section 2.4 . The pairwise re duction is performed based on the best matches of identified cor responding regions between the fragments by analysing the simi larities of curvatures in scalespace ( Section 2.5 ). The merged frac ture lines after the pairwise reduction is calculated as described in Section 2.6 . The scalespace curvature matching is repeated un til all fragments are reduced. The set of fracture reconstructions produced from the second module ( Section 2.7 ), are used as input for the third module of the algorithm. Finally, the algorithm deter mines the best solution based on all performed reconstructions as described in 2.7.1.

2.2. Generation of triangular surface models

The proximal humeral anatomy is composed of an (outer) corti cal layer and (internal) cancellous bone. The cortical layer is com pact and dense and appears bright in a CT image. The cancellous bone is a porous structure, which is less dense than the cortical layer. Our reconstruction algorithm is based on surface models that represent only the cortical layer and not the cancellous bone of the humeral fragments. For the clinical cases, the generation of the triangular surface models was performed by an experienced orthopaedic surgeon with clinically applied segmentation methods ( Murase et al., 2008; Vlachopoulos et al., 2016b ). Thereby, we used thresholding, man ually correction of connected fragments, region growing, and the marching cubes algorithm ( Lorensen and Cline, 1987 ). For the cadaver experiments we used the triangular surface models provided by the authors of the study of Fürnstahl et al. (2012) which used the bone enhancement fil ter of Descoteaux et al. (2006) for the segmentation of the cortical layer. Best results were achieved with a filter range of 3.5 mm, which was adapted to match the maximal cortical thickness of the proximal humerus. The average cortical thickness towards the glenohumeral joint was evaluated to be between 0.75 mm and 1 mm.

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 145

Fig. 2. Overview of the reconstruction of a proximal humeral fracture.

Fig. 3. Fracture surface representation. The fracture surfaces of the tuberosity fragment and the humeral shaft fragment are simplified by a fracture line, represented by a sequence of points (denoted by red spheres). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.3. Fracture surface representation

The thickness of the cortical layer is about 4 mm in the proxi mal shaft region and decreases towards the joints ( Fürnstahl, 2010; Skedros et al., 2016 ) As a consequence, the fracture surfaces are extremely narrow, and can be can be considered as a narrow path. Therefore, our idea was to represent the fracture surfaces by a frac ture line along the path. Formally, a fracture surface FS of a given fragment is simplified by a fracture line FL, i.e., a sequence of points p i ∈ FS . The anno tation of the points p i ^{on the triangular surface model was per} formed manually by setting densely sampled points as illustrated by spheres in Fig. 3 . The points were projected on the fracture sur face and centred between the inner and outer contour of the cor tical layer. To ensure equally distanced points, the algorithm per formed a cubic spline interpolation. The distance between the in terpolated points was selected to be r F L ^{= 1 mm. The order of the} point sequence specifies the direction of the path. Therefore for each fragment, the fracture surface was represented by a sequence of points in clockwise and counterclock wise direction.

2.4. Curvature calculation

The calculation of the curvature of a 3D curve is well known from differential geometry ( Kreyszig, 1959; Salomon, 2007 ). How ever, it is also known that curvature computation from noisy data is problematic which makes the application to the reconstruction

of fractured bones very difficult. To tackle this problem we in troduce a scalespace representation of the curvature of the frac ture lines as follows. We define the scaledependent local shape of FL around a point p i ^{by an approximating cubic Bspline with} one polynomial piece B [ −d,d ] p i . The scale parameter d ≥2 denotes the length of the fracture line segment { p s } ∈ FL around p i ^{to be consid} ered, where s = { i −d, . . . , i, . . . , i + d } . The Bspline is constructed in a weighted meansquare sense as described in De Boor (1978 ), minimizing



s w s ^(B [ −d,d ] p i ( p s ) −p s ) 2 , (1)

where B [ −d,d ] p i ( p s ) is the value of the spline B [ −d,d ] p i for point p s ^and the weights are selected as

w s ⁼

 j ² , s < i j 4 , s = i ( 2 d + 2 −j ) 2 , s > i and j = { 1 , 2 , . . | s | }

The weighting function enforces that the approximating cubic Bspline calculated around a point p i ^{passes through the spline’s} centre point p i . ( Fig. 4 a). The set of approximating cubic Bsplines calculated for scale d represents the fracture line FL of the frag ment in scale d ( Fig. 4 b). In other words, the entire fragment is characterized by the fracture line.

146 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Fig. 4. Curvature Calculation. a) Approximating cubic BSpline are illustrated for scales d = 5 (black), d = 10 (blue), and d = 15 (green) centred on a point p i ^{of the fracture} line. b) The Bsplines for all p i for scale d = 10 represent the fracture line of a fragment. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The curvature vector k B [ −d, d ] p i of B [ −d,d ] p i is defined analog to the

standard definition of the curvature vector of a parametric curve as described in De Boor (1978 ). The local shape of a fracture line can be classified into concave or convex parts, depending on the orientation of the curvature. Similar as for a jigsaw puzzle, corresponding border segments of two fragments have opposing curvature orientations. To take these characteristics into account for the matching of corresponding cur vatures, we introduce the signed curvature value around a point p i

k s B [ −d, d ] p i = ∥ k B [ −d, d ] p i ^{∥ sign}  k B [ −d, d ] p i · v p i  . (2)

where · denotes the dot product and v p i ^{the surface normal of the} underlying triangle mesh at p i . For two signed curvature values ks 1 ^{and ks} 2 ^{we define the cost} c f( k s 1 , k s 2 ) = k s 1 ^{+ k} s 2 , i.e., for two ideal matching fracture frag ments the cost function is zero for all curvature values, since the absolute curvature values are equal but their signs are opposite. The cost function will be used for the selection of matching candi dates as described in the following Section 2.5 .

2.4.1. Scale space representation The scale space representation is used to search for matching candidates between two fracture lines in the same scale. For a scale d the curvature along the fracture line can be expressed by the curvature function

K F L ( p i, d ) = k s B [ −d, d ] p i , ∀ p i ^{∈ F} L. (3)

The scalespace D is generated for all fragments by calculating the curvature function of Eq. (3) for all p i ^{and scales d} ∈ {5, 10, 15, 20, 25}. The scalespace representation of the humeral shaft frag ment is illustrated in Fig. 5 .

2.5. Pairwise reduction

Given two different fragments F l ^{and F} m ^{we seek for the best} possible pairwise reduction(s). We perform the selection of the best pairwise matches in two steps. In the first step we deter mine corresponding regions of the curvature of the fracture lines FL l and FL m ^{where the cost function cf is minimal. The correspond} ing regions are detected in the same scale of the scalescape D l ^and D m ^{as described in 2.5.1} . Based on the corresponding regions, we

Fig. 5. Scalespace representation of the curvature function. The sign of the curva ture value is encoded by the colour (concave: blue, convex: red) and the periodicity is depicted as dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

perform the pairwise reductions and evaluate the reductions which result to the smallest reduction error as described in 2.5.2 . For each combination of two fragments we will keep the best possible pairwise reductions for the fracture reconstruction as described in 2.7 .

2.5.1. Corresponding regions For two fracture lines the curvature functions of Eq. (3) are analysed to identify corresponding regions in each scale. Only highly promising candidates have to be investigated to reduce the number of possible combinations of matching curvature regions. The clinical observation indicated that especially the fragment of the greater tuberosity has regularly a distal triangular tip, which

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 147

Fig. 6. Visualization of matching intervals. The interval is illustrated on scale d = 20 of the curvature functions of FL l ^{and FL} m ^{with corresponding peaks pks} u ^{and pks} v and lengths λ1 and λ2 .

matches to a corresponding defect in the shaft ( Gerber et al., 2004 ). In the surgery, the typical first step is to reduce the greater tuberosity ( Gerber et al., 2004 ). These findings led us to the as sumption that peaks in the curvature function represent promising candidates, as a prominently shaped region of the fracture surface is supposed to have high curvature. Therefore, our approach is to analyse the peaks pks ( FL , d ) of each curvature function K FL . Peaks are detected by calculating local maxima of the absolute curvature function | K FL ( p i , d )|. To further reduce the number of candidates, all local maxima that are smaller than 10% of maximum curva ture value (i.e., max | K FL ( p i , d ) | ∀ p i ∈ FL ) can be rejected. Thereafter, in each scale corresponding pairs [ pks u , pks v ] of peak points are determined for ∀ pks u ∈ pks ( FL l , d ) and ∀ pks v ∈ pks ( FL m , d ), where K F L l ( pk s u , d ) K F L m ( pk s v , d ) < 0 . Candidate selection along the fracture lines FL l ^{and FL} m can then be performed for each scale d and corresponding pairs [ pks u , pks v ] by finding the largest intervals { p q . . . p q + λ^} and { p r . . . p r+ λ^{} of length λ, with pk} s u ^{∈ {} p q . . . p q + λ^} , p k s v ^∈ { p r . . . p r+ λ^{} and λ ≥10, for which it holds that}

 c f  K F L l ( p q + t , d ) , K F L m ( p r+ t , d )   < ^| K F L l ( pk s u , d ) | + | K F L m ( pk s v , d ) |

2 ∀ 0 ≤t ≤λ .

(4)

It should be mentioned that whenever p q + λ ^{or p} r+ λ ^{exceed the} domain of K FL , the function K FL ^{is assumed to be periodic and, in} addition, it holds that λ < min (| FS l |, | FS m | ) where | | denotes cardi nality. According to Eq. (4) matching intervals can be expressed by the peak points pks u ^{and pks} v ^{if two lengths λ}1 and λ2 are introduced. For example, in Fig. 6 , the matching inter val ci = [ pk s u , pk s v , λ1 , λ2 ] between two fracture lines equals the set of points { pk s u −λ1 ^. . . p k s u . . . p k s u + λ2 ^{} = {} p q . . . p q + λ^{} and} { pk s v −λ1 ^. . . p k s v . . . p k s v + λ2 ^{} = {} p r . . . p r+ λ^} . The matching inter vals will be used to calculate the relative transformations, that align the fragments F l ^{and F} m ⁽ Section 2.5.2 ), and in addition to evaluate the pairwise reduction.

2.5.2. Evaluation of the pairwise reduction The method of calculating the pairwise reduction was inspired by the surgical technique for the alignment of two fragments. Here, after identification of matching features ( Gerber et al., 2004 ), the fragments are manipulated by the surgeon to align the fragments as accurately as possible. The pairwise reduction is achieved by the calculation of the transformation which aligns the two fragments. In most of the previous approaches ( JimenezDelgado et al., 2016 ), which use a ICPbased algorithm to perform the reduction, the correspondence between the surface points is not known and has to be deter

mined by iterative calculation of correspondences. In contrast to simple ICPbased approaches, the correct correspondence between points is implicitly given by our method ( Fig. 7 ). The rigid trans formation T ci ^{which aligns the point sets {} p pk s v −λ1 ^. . . p pk s v + λ2 } and { p pk s u −λ1 ^. . . p pk s u + λ2 } , is calculated using absolute orienta tion ( Horn, 1987 ). However, the 3D orientation of fracture lines necessitates introducing a second measurement. Previous studies from archaeology used the torsion as a signature of the 3Dcurve ( Papaioannou and Theoharis, 2003; Üçoluk and Toroslu, 1999 ). Our approach was to evaluate the distance between the points of the fracture lines after the reduction. Papaioannou and Theo haris (2003) considered two fragments surfaces as matching can didates if the length of the boundary segments was at least one fourth of the arc length of the shortest boundary. As we aimed to avoid a fixed threshold value, our approach was to assign a higher weight to larger matching intervals. Therefore, we defined the pair wise reduction error measure prem , by the root mean squared er ror (RMSE) of the Euclidean distance between the aligned point sets, normalized by their arc length:

pre m ci ^{= RMSE}

λ (5)

Assuming that two matching intervals of different sizes would result in the same RMSE , the normalization by the arc length would favour the larger intervals. Hence, the set of matching intervals for the fracture lines FL l and FL m ^{(as calculated in Section 2.5.1} ) is

C I F L l , F L m = { ci | pk s u ^{∈ pks} ( F L l , d ) , pk s v ^{∈ pks} ( F L m , d ) , ∀ d ∈ D }

For each matching interval we calculate the rigid transforma tion that aligns the corresponding point sets of the fracture lines. Correspondingly the calculated set of transformations to reduce the fragments F l ^{and F} m ^{is defined as}

T F L l , F L m ⁼ T ci | ci ∈ C I F L l , F L m

For the further iterative fracture reconstruction, we merge the fragments F l ^{and F} m ^{and calculate merged fracture line(s) as de} scribed in 2.6 .

2.6. Calculation of the merged fracture line

For the determined set of transformations T F L l ,F L m ^{we calculate} merged fracture lines used for the further fracture reconstruction as illustrated in Fig. 7 . Without loss of generality it holds that |{ p i }| > |{ q i }|, p i ∈ FL l , q i ∈ FL m ⁽ Fig. 7 , a). The steps of the algorithm for the calculation of a merged frac ture line are as follows:

1. Apply the transformation T ci ^{to the sequence of points of the} fracture line FL m ⁽ Fig. 7 , b). 2. Calculate the residual points of the fracture lines F L res _ l ⁼ F L l ^{\ {} p pk s u −λ1 ^. . . p pk s u + λ2 } and F L res _ m ⁼ F L m ^\ { p pk s v −λ1 ^. . . p pk s v + λ2 } ( Fig. 7 , c). 3. Determine the end points ( p start _ l ^{and p} start _ m ) and ( p end _ l ^and p end _ m ) of F L res _ l ^{and F} L res _ m ^{on both sides.} 4. Connect F L res _ l ^{and F} L res _ m ^{by cubicspline interpolation be} tween p end _ l ^{and p} end _ m ^{and between p} start _ l ^{and p} start _ m ( Fig. 7 , d), and subsample the connection with a sampling size of r FL ^{by introducing additional points {} p co n start } and { p co n end } . 5. Determine the point p merge ^{of F} L res _ l ^{with the maximal distance} to the nearest point of F L res _ m ⁽ Fig. 7 , e). 6. Starting from p merge ^{create a directed weighted graph of} { F L res _ l ^{∪ p} co n start ^∪ F L re s m ^{∪ p} con end } , where each point is connected to all subsequent points within a radius of r graph ^{= 1} . 5 mm and the edge weights correspond to the distance between the points of an edge. In order to obtain a similar distance r FL , as defined in Section 2.3 , between the

148 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Fig. 7. Calculation of the merged fracture line. a) The calculated transformation T ci , which aligns the point sets of the matching interval (denoted by cyan spheres) is used to b) align FL l ^{and FL} m . c) The point sets of the matching interval are removed. d) The residual points of the fracture lines F L res _ l ^{and F} L res _ m ^{are connected with p} co n end ^by cubicspline interpolation between p end _ l ^{and p} end _ m . e) Starting from p merge ^{create a directed graph and use the shortest path algorithm to determine the merged fracture} line FL l, m . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Importance of the order of pairwise reductions, illustrated for a proximal humeral fracture with five fragments. 1) The shaft fragment (beige) is first merged with the blue fragment and 2) the cyan fragment is merged with the orange fragment. 3) Thereafter, the results of 1) and 2) are merged and 4) finally, the result of 3) is merged with the head fragment (yellow). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

sequence of points of FL l, m , r graph ^{has to be greater than r} FL and smaller than 2 ∗r FL . Therefore, we set r graph ^{= 0} . 5 ∗( r F L ⁺ 2 ∗r F L ) . 7. Calculate the merged fracture line FL l, m ^{using the shortest path} algorithm ( Dijkstra, 1959 ) ( Fig. 7 , f).

2.7. Fracture reconstruction

Given a fracture with n fragments { F 1 , F 2 , . . . , F n } the goal of the fracture reconstruction algorithm is to calculate n correspond ing transformations, which will reduce the fragments to their pre traumatic position. The fracture reconstruction is performed by evaluating the reduction pairwise and repeating the process iter atively for merged fragments as described in 2.3 to 2.6 . The order of pairwise reduction is crucial for the outcome of any fracture re construction algorithm. For example, the correct reduction of one fragment might only be possible if the pairwise reduction of two other fragments has already taken place. Fig. 8 demonstrates the relevance of the reduction sequence for a proximal humeral frac ture with five fragments. Because the optimal order cannot be determined a priori, our algorithm processes all possible combinations and determines fi nally the best solution by evaluation of the reconstruction re sults as described in 2.7.1. For a threepart fracture, it is neces sary to calculate three combinations of pairwise reduction, which corresponds to 3 x 2 = 6 pairwise reductions in total. For a four

part fracture, there are fifteen combinations, which yield in total 15 x 3 = 45 pairwise reductions. For a fivepart fracture, there are 135 combinations, which yield in total 135 x 4 = 540 pairwise reduc tions. Fig. 9 illustrates the sequence of combinations. For the fracture reconstruction we consider only the most promising b max ^{= 20 pairwise reductions, which yield an ac} ceptable combination effort. To ensure that we initially consider matching intervals in each of the scales, the b max ^{pairwise reduc} tions are determined as follows. For each scale we calculate the b max = 20 best matching intervals with the smallest prem (5) and keep, thereafter, the overall b max = 20 best matching intervals with the smallest prem . Thereby, the fracture reconstruction is performed in iterations. In the first iteration, the b max ^{promising pairwise reductions of two} fragments are determined. The curvature and the scale space of the merge fracture lines are calculated. In the next iteration, the b max ^{pairwise reductions between the merged fragments and the} next fragment of the reduction sequence are determined. The over all maximum number of calculations for each combination is given by N max _ per _ comb ^{= b} max n −1 . After the last iteration, the algorithm determines the best solution as described in 2.7.1 below.

2.7.1. Best solution Lastly, the algorithm determines the best solution based on all performed reconstructions. The calculation of the total fracture

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 149

Fig. 9. Combination of pairwise reduction for the three and fourpart fractures.

reconstruction error measure frem is based on the finally aligned points of the fracture lines of all fragments after the last iteration. In contrast to the evaluation of the pairwise reductions, the corre spondences between the point pairs of the fracture lines are not known. We calculate frem as:

frem = 1  n l=1 ^| F l l |

n 

l=1

| F L l | 

i =1 min ( f dist ( p i , Q ) ) , p i ^{∈ F} L l ,

Q =

q = n 

q =1 , q ̸ = l F l q (6)

where f dist ( p i , Q ) returns the Euclidean distances between a point p i ^{and all points of the point set Q} . The fracture reconstruction, which yielded the smallest frem value, was regarded as the best solution. We have to note that the automatic selection of the best solu tion might fail if only a partial reconstruction is performed, i.e., if relevant fragments of the proximal humerus are not considered in the reconstruction algorithm. One possible reason for this could be that the surgeon failed to identify and segment all fragments. For example, a fracture reconstruction which leads to a folding of the fragments could have in these cases a smaller frem value than the optimal reconstruction. A warning mechanism was implemented in the automatic se lection for the best solution to inform the user that relevant in formation might be missing, that are required to achieve the re construction. The idea of the mechanism is to verify whether the surface area area fragments ^{of all considered fragments is comparable} to the surface area of a healthy proximal humerus area expected . The mechanism was only necessary if a humeral head fragment was present. The sum of the surface area of all fragments area fragments was calculated as described in Section 3.2 . The surface area of the

shaft fragment was excluded since the size depends mainly on the acquisition length as described in 3.2. The approximated area area expected ^{was estimated by the surface area of a sphere fitted to} the surface points of the humeral head fragment (see Fig. 10 ). In case of a fracture with a humeral head fragment the algorithm calculates the expected surface area area expected ^and area fragments . If area expected ≥2 ∗area fragments ^{the algorithm raises a} warning, that the best solution might not be the solution with the smallest frem . In these cases the solutions have to be inspected by the user.

3. Results

3.1. Datasets

A consecutive series of eight patients, with a proximal humeral fracture treated with open reduction and internal fixation (ORIF), in the department of orthopaedics at the University Hospital Bal grist, Zurich, Switzerland, between January 2016 and July 2016, were used for the clinical evaluation. The CT scans were obtained according to a standard scanning protocol used for the preop erative evaluation of proximal humeral fractures (slice thickness 1 mm; 120 kV; Philips Brilliance 40 CT, Philips Healthcare, The Netherlands). The cantonal ethics committee of Zurich, Switzerland approved the study (KEKZHNr. 20130586). To compare our ap proach with the stateoftheart method ( Fürnstahl et al., 2012 ), we applied our reconstruction algorithm also to data previously pub lished by Fürnstahl et al. (2012) , i.e., four clinical cases and four ar tificially created fractures on cadaveric humeri. In total we applied our method to sixteen complex fractures of the proximal humerus.

150 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Fig. 10. Approximation of the expected area of reconstruction. The surface area of a sphere, fitted to the surface points of the humeral head fragment, was used to approximate the expected surface area of reconstruction.

3.2. Fragment characteristics

During the reconstruction of proximal humeral fractures by ORIF the larger fragments are reduced, while smaller fragments are left in place as local bone grafts to facilitate the healing pro cess. Therefore, previous reconstruction algorithms automatically removed small fragments before the reconstruction. One measure of identifying small fragments is the area of the mesh triangles, i.e., all fragments with a summed area of the triangles smaller than 10 0 0 mm ² were small and not considered in the reconstruction in the study of Fürnstahl et al. (2012) , since contralateral matching did not robustly work for them. In the current study, we first analysed the size of all fragments of the clinical cases that were classified by an orthopaedic surgeon and considered as being large enough for ORIF. To measure the size of a fragment, we calculated the oriented bounding box (OBB) from all model points. The OBB was previously described as a valuable method for the 3D measurement of the humeral length for clin ical applications ( Vlachopoulos et al., 2016a ). The area of a frag ment was approximated by the product of the two longest sides of the OBB. The characteristics of all thirty fragments, that were classified by the surgeon as large fragments and used for ORIF, are summarized in Table 1 . The shaft fragments are not included, as the size is arbitrary and mainly depend on the scanned length of the upper arm during CT acquisition. The area of the smallest fragment of these thirty fragments was 195 mm ². Therefore, we used for the fracture reconstruction of the clinical cases and of the cadaver experiments all fragments with an area greater than are a min ^{= 195 mm} ².

3.3. Clinical evaluation

All reconstructions of the twelve clinical cases resulted in an ac curate approximation of the pretraumatic anatomy. The quality of the best solution provided by the algorithm was evaluated by two

Table 1 Characteristics of the fragments of the clinical cases. The area of the fragments is summarized, which were con sidered by the surgeon as relevant during fracture recon struction. The intermediate fragments are all fragments excluding the head fragments.

Area (mm ²)

Mean SD Range

Head (7) 2320 682 1812 −3319 Intermediate (23) 880 520 195 −1954 All fragments (30) 1216 828 195 −3319

Table 2 Characteristics of the fragments of the cadaver experiments. The area of all fragments is summarized, which were used for the fracture reconstruction.

Area (mm ²)

Mean SD Range

Head (4) 2925 228 2749 −3232 Intermediate (10) 860 486 364 −1958 All fragments (14) 1409 1037 364 −3232

surgeons, trained in orthopaedic surgery and specialized in shoul der surgery. The evaluation was based on a presentation showing the 3D reconstructed fracture from multiple, standardized view points. The examiners were told to evaluate the proposed recon struction, i.e., whether the reconstruction would restore the pre traumatic anatomy or whether a malposition of the fragments was apparent. In clinical studies the postoperative malposition is nor mally assessed on conventional radiographs and a malposition of a fragment is assumed, when a displacement of at least 1 cm or of 45 ° is present ( Gerber et al., 2004 ). We defined for the clinical evaluation a much lower cutoff of 5 mm or 10 °, respectively. For the evaluation, we used the following 5 pointscale:

1. Acceptable without modification 2. Acceptable after small modifications 3. (i.e. correction of one fragment < 10 ° or < 5 mm) 4. Acceptable after large modifications 5. (i.e. correction of more than one fragment or one fragment > 10 ° or > 5 mm) 6. Not acceptable at all 7. Determination of the correct alignment is not possible

In Fig. 11 , we present the reconstruction results of all clini cal cases sorted by the number of fragments. Both surgeons rated eleven of the twelve clinical cases as acceptable without modifica tion. Case 6 was rated by both surgeons as acceptable after small modifications, since the reconstruction resulted in a small malro tation of the humeral head.

3.4. Quantitative evaluation based on cadaver experiments

The accuracy of our method was evaluated based on four cadaveric data sets, previously acquired in the study of Fürnstahl et al. (2012) . The CT scans had the same axial resolution of 1 mm, as the clinical cases. The fragment size of the cadaver bone fractures was also comparable ( Table 2 ). All cases could be fully reconstructed. The best solution provided by the algorithm is illustrated in Fig. 11 . In case 13, the warning mechanism of the algorithm indicated that the proposed solution with the smallest frem error might be not the best one. The reason was that the creation of the artificial fracture had resulted in a comminuted fracture with a large number of small fragments. These fragments were not included into the fracture reconstruction algorithm, simi

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 151

Fig. 11. Reconstruction results. The cases are sorted by the number of fragments. Clinical Cases: Case 1 (two fragment), Case 2–6 (three fragments), Case 7–11 (four fragments) and Case 12 (five fragments) Cadaver Experiments: Case 13 (three frag ments) and Case 14–16 (five fragments).

lar as in the previous study ( Fürnstahl et al., 2012 ). In this case the area fragments ^{was 0.46} ∗area expected . The goal of a fracture reconstruction algorithm is to reduce the fragments such that the pretraumatic anatomy is approximated as well as possible. The approximation of the pretraumatic anatomy can be quantitatively measured in the cadaver experiments by comparing the difference of the reconstructed surface of the frag ments with the surface of the original unfractured bone model. Here, the distance error was defined as the Euclidean distance between all points of the outer corticalis of the reduced fragments to the closest point on the outer corticalis of the original bone model ( Fig. 12 ). The distance errors of all cases reconstructed with the presented algorithm and the reference method ( Fürnstahl et al., 2012 ) are presented in Table 3 .

3.5. Order of pairwise reduction

The algorithm found the best solution in fourteen of the eigh teen cases by subsequent merging of the fragments (similar as il lustrated in the top rows of Fig. 9 for the fourpart fractures). In four cases, the algorithm found the best solution by performing first two pairwise reductions of two fragments (i.e., F 12 ^{and F} 34 ), before the final reduction was completed (similar as illustrated in the bottom rows of Fig. 9 for the fourpart fractures and in Fig. 8 ).

Fig. 12. Distance error. Euclidean Distances are calculated for all surface points of the outer cortical layer between the reduced fragments and the original bone model.

3.6. Runtime

All parts of the proposed algorithm were implemented in

MATLAB (2015) . The runtime scales with the number of fragments of the fracture. Table 4 summarizes the runtimes of our algo rithm (curvature calculation, pairwise reduction, calculation of the merged fracture line and fracture reconstruction). To compare the results with the reference method, we summarized the runtime of the reference method for the contralateral assembly and the fracture surface assembly, without the preprocessing time (fracture segmentation and fracture surface extraction).

3.7. Comparison with statisticalshapemodel based fracture reconstruction

We compared the proposed method with a StatisticalShape Model (SSM) based fracture reconstruction approach, which may be considered as an alternative method. For the prediction of the proximal humeral anatomy we used the whole distal humeral model for the fitting of the statistical shape model. As standard scanning protocols used for the preoperative evaluation of proxi mal humeral fractures do neither include the whole humeral shaft nor the distal part of the humerus, we could not apply the SSM based fracture reconstruction to the original data of Cases 1 to 16. Therefore, we performed experiments on simulated simple frac tures of humeral models obtained from healthy cadavers as illus trated in Fig. 13 . We used 3D triangular meshes of 50 right cadav eric humeri without a pathological condition. The segmentation of the humeri was performed in an automatic fashion using a previ ously described segmentation algorithm ( Gass et al., 2014 ). The 3D triangular meshes of the 50 right humeri were brought into cor respondence using a nonrigid registration algorithm ( Lüthi et al., 2016 ). Subsequently, all meshes were rigidly aligned using Pro crustes alignment ( Umeyama, 1991 ) to one humeral model. A SSM of the aligned meshes was built by performing a Principal Compo nent Analysis (PCA) ( Jolliffe, 2002 ). We simulated a simple fracture with one fragment of the prox imal humerus. The length of the humeral head fragment was de fined to be 15% of the length of the whole humerus. Subsequently we fitted the SSM to the distal 85% of the humerus to predict the proximal humeral anatomy. In total, we analysed the simu lated fractures on 15 cadaveric humeri with leaveoneout cross validation. We refer to ( Albrecht et al., 2013 ) for the mathematical details of the prediction procedure. The distance error between the original bone model and the predicted bone model was calculated as for our proposed method in the experiment before ( Table 5 ). For the SSM experiment, we assume a simple fracture with one frag ment that was also not present in the cadaver experiments.

152 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Table 3 Comparison of the distance errors between the proposed algorithm and the reference method of Fürnstahl et al. (2012) . Mean, standard deviation and range of the error are given.

Cadaver Distance error (mm)

Case Fragment Presented method Reference method

Mean SD Range Mean SD Range

13 Head 1.31 0.73 0.02–4.99 1.81 0.99 0.02–6.02 Intermediate 1.09 0.57 0.02–2.87 1.70 0.83 0.08–3.79 14 Head 1.68 0.89 0.02–4.55 0.89 0.38 0.04–4.32 Intermediate 0.76 0.41 0.02–2.75 0.90 0.34 0.04–2.70 15 Head 1.18 0.56 0.03–4.69 0.58 0.30 0.02–3.87 Intermediate 0.95 0.44 0.03–2.84 1.16 0.56 0.01–3.68 16 Head 1.36 0.74 0.03–3.94 1.48 0.77 0.02–3.45 Intermediate 0.95 0.67 0.02–3.88 0.80 0.48 0.02–3.84

All Head 1.41 0.78 0.02–4.99 1.23 0.85 0.02–6.02 Intermediate 0.90 0.56 0.02–3.88 0.97 0.56 0.01–3.84

Fig. 13. Statistical Shape Model. a) Humeral models without a pathological condition were used for the simulation of a simple proximal humeral fracture with one fragment. b) Based on the distal 85% of the humerus, we predicted the proximal humeral anatomy with the SSM. c) d) The distance error between the original bone model and the predicted bone model was calculated.

Table 4 Average runtime of the presented method implemented in MATLAB in sec onds (3.4GHz Intel Core i72600 CPU, 16GB RAM) compared to the runtime of the CPU version of the reference method (P4 3.2GHz).

Runtime (s)

Number of fragments Presented method Reference method (CPU)

3 121 2859 4 4220 5742 5 11,322 9127

Table 5 Distance error for the simulated simple fractures of the SSMbased method. Mean, standard deviation and range of the error are given.

Distance error (mm)

Mean SD Range

Cadaver ( n = 15) 1.55 0.74 0.05 −7.79

4. Discussion

In this work, we presented a novel method for the fully au tomated anatomical reconstruction of proximal humeral fractures.

One major advantage of the presented method in contrast to the stateoftheart ( Fürnstahl et al., 2012 ) is that the knowledge of the contralateral anatomy is not necessary anymore. Therefore, our approach is even applicable in presence of bilateral pathological conditions. At the same time, the acquisition of the contralateral anatomy necessitates an additional CT scan, leading to increased radiation exposure for the patient. This is a major hurdle for the clinical introduction of a new method in the hospital. Our algo rithm tackles this problem by performing the fracture reconstruc tion solely based on the information of the fractured surfaces with out relying on a reconstruction template. Certainly the fractured surface of a fragment does not represent the whole morphology of the fragment. However, we demonstrated on twelve clinical cases that the information of the fractured surfaces is already sufficient to accurately reconstruct proximal humeral fractures with respect to clinical requirements. Even if it is possible to replace the contralateral anatomy by a different reconstruction template, e. g., a statistical shape model as in Albrecht and Vetter (2012) , the fragments at the proximal humerus are often considerably displaced and malrotated, making their registration to the template difficult ( Fürnstahl et al., 2012 ). According to our experience, a reconstruction using a SSM can not be reliably performed if data of the distal part of the humerus

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 153

Fig. 14. Influence of a plastic deformation of a fragment on the ICP registration result. a) The fragment shows a plastic deformation and the distance error after ICP registration b) The deformed area of the fragment was removed before the ICP registration.

is missing. Unfortunately, in almost all clinical CT scanning pro tocols for proximal humeral fractures neither the whole humeral shaft nor the distal humerus is included. The reason for omitting a scan of the distal humerus is the reduction of radiation exposure. It is a clear advantage of our method that radiation exposure is kept minimal. Furthermore, Albrecht and Vetter (2012) pointed out than the length of the bone has to be correctly initialized in order to achieve an accurate fracture reconstruction, which further limits the application of their method. We demonstrated, that even if a simple fracture with one fragment is present, the distance error of the predicted model is higher that with our proposed method. Most of the current approaches for fracture reconstruction have in common that an ICPbased algorithm is used to perform the re duction of the fragments ( JimenezDelgado et al., 2016 ). The ten dency of the ICP to fall into local minima might be particularly problematic in case of the proximal humerus due to its almost spherical shape ( Fürnstahl et al., 2012 ). Furthermore, an ICPbased algorithm might introduce an error in the reconstruction of the pretraumatic anatomy in presence of a plastic deformation of the fragments. The relevance of threedimensional analysis for the di agnosis of plastic deformation of a presumed intact intraarticular surface of a complex malunited humeral fracture has already been emphasized ( Vlachopoulos et al., 2016b ). The influence of the plas tic deformation onto the ICPregistration results is demonstrated in Fig. 14 . Although only a small part of the fragment is deformed, the ICPregistration yields an incorrect result if a plastic deforma tion is present. In the first column of Fig. 14 the tuberosity frag ment of Case 14 is presented as it was used for the fracture recon struction and after registration onto the unfractured humerus. In the second column, the plastic deformed part of the fragment was first removed before the ICPregistration. It can be clearly seen how the inherent local deviation caused by the plastically deformed re gion deteriorates the reconstruction quality of the entire fragment. On the contrary, the fracture reconstruction with our presented method is not compromised as it did not rely on the fracture line segment affected by the plastic deformation.

Compared to other fracture reconstruction approaches, that per form the reconstruction based on the fracture surfaces of the fragments our algorithm has clear advantages but also some drawbacks. The main drawback of our approach compared to the method of Okada et al. (2009) or Kronman and Joskow icz (2013) is the significantly higher runtime. Nevertheless, these registration methods can only handle fragment displacements up to a certain degree. Okada et al. (2009) limited the displace ments to 30 ° and 20 mm and Kronman and Joskowicz (2013) lim ited the displacement to 20 ° and 30 mm around each axis. The poseinvariant feature of our algorithm is therefore of particu lar importance for fractures of the proximal humerus, due to the observed greater displacements and a great advantage com pared to the previous algorithms. Furthermore, the method of Okada et al. (2009) could register the fracture surfaces only if the order of alignment as well as the corresponding regions were cho sen beforehand. Best result were obtained with a method that combined fracture surfaces and the contralateral anatomy, and Okada et al. (2009) pointed out, that constraints posed by frac ture surfaces alone were insufficient to perform registration with acceptable accuracy. Okada et al (2009) and Zhou et al. (2009) de termined the order of pairwise reduction in beforehand. However, such a strategy may not be possible for all fractures of the proxi mal humerus. Already in our dataset we observed four cases where the strategy of predefining the order of matching by the fragment size would not work (illustrated in Fig. 8 ). An additional limita tion of Kronman and Joskowicz (2013) is that the method is only applicable for bone fractures consisting of two fragments and the method was not extended to fractures with multiple fragments. The quantitative evaluation based on cadaver experiments demonstrated that the reconstruction results of our method are almost equal for head fragments and better for the small fragments, where the reconstruction is even more challeng ing, compared to the stateoftheart ( Fürnstahl et al., 2012 ). The mean distance error of reconstructed intermediate frag ments was 0.90 mm ± 0.56 mm for our method compared to

154 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

0.96 mm ± 0.56 mm for the stateoftheart. The distance error of the head fragments was 1.41 mm ± 0.78 mm compared to 1.23 mm ± 0.85 mm for the stateoftheart. As the fracture reconstruction is based on the fracture lines, the reconstruction of comminuted fracture might be more prob lematic in some cases. Although the fracture pattern, particularly of the cadaver experiments, showed a variable amount of fracture comminution, the fracture reconstruction was still possible in the presented cases. However, a comminuted fracture that would al ter more the corresponding fracture lines of two fragments might be more difficult to handle with the presented method compared to an algorithm that relies on a reconstruction template, similar as the method of Fürnstahl et al. (2012) . One elegant solution in cases of comminuted fractures would be to include user interaction in the optimization process. For instance, a haptic surgical planning system ( Fornaro et al., 2010; Harders et al., 2007; Kovler et al., 2015; Olsson et al., 2013 ) could be used to make small modifica tions to automatically reduced fragments based on additional clin ical considerations. The computational efficiency of a reconstruction algorithm is a clinically important factor, as it significantly contributes to the overall preparation time for the surgery. For a fracture with four fragments, the runtime of our method was 4220 s, which outper forms the CPUversion of the reference method (5742 s). For a frac ture with three fragments the runtime benefit was much greater (121 s compared to 2859 s) when our algorithm is used. Fractures up to three or four parts are treated by ORIF, but a surgical recon struction of fractures with more parts is very unlikely to be suc cessful ( Gerber et al., 2004; Wijgman et al., 2002 ). There would be several options to reduce the complexity of the calculation and therefore the runtime. The implementation of the algorithm in a compiled language, or even better on a graphical processing unit, would allow reducing the runtime (i.e., the run time of the algorithm of Fürnstahl et al. (2012) was reduced about a factor of 20). In the presented algorithm, the order of the point sequence specifies the direction of the path of the fracture lines. Therefore, for each fragment, the fracture surface was represented by two sequences of points, one in clockwise and one in counter clockwise direction. We performed the pairwise matching calcu lated from the sequence of points of one fragment in one direction and for the other fragment in both directions. The automatic de tection of the outer contour of the cortical layer would permit to determine the appropriate direction, and would halve the calcula tion of pairwise matching. Additional speedup would be possible, if the number of matching candidates per pairwise reduction can be reduced already in the early stage of the algorithm, i.e., by in corporation of geometric constraints of the merged fragments into the decision. Our method relies on some heuristically defined parameter val ues which we shortly discuss. For the scalespace representation the scales d ∈ {5, 10, 15, 20, 25} were considered. The lower range was set to scale 5, as significant shape features on the fracture bor ders had a length of at least 10 mm length. Thereby, a subset of at least 10 consecutive points influenced the calculation of the curva ture at each point p i . The minimal area of a fragment, which was considered by surgeons as relevant for the fracture reconstruction, was 195 mm 2 . The arclength of such a fragment would be approx imately 50 mm. The upper value of d was set to d = 25 to consider all 50 points of the smallest fragments in the calculation of the cur vature. The maximal number of best matches, which were selected after evaluation of the pairwise reductions for the further calcula tions in each iteration was set to be 20. The number was defined to be relatively high in order to enable an automatic fracture re construction for all cases. The presented approach assumes that the fracture surfaces can be represented by fracture lines. The algorithm might, therefore,

be applicable for further anatomies, if this representation is also possible. Okada et al. (2009) also uses fracture lines as an rep resentation and, therefore, we strongly believe that our algorithm can be extended to other anatomies such as diaphyseal and meta physeal fractures of long bones, i.e. femoral fractures as presented by Albrecht and Vetter (2012), Kronman and Joskowicz (2013) or Okada et al. (2009) .

5. Conclusion and future work

In this paper, we proposed an algorithm for the reconstruction of complex proximal humeral fractures. The key idea of the ap proach is the use of the curvature scalespace for matching char acteristic features between the fragments. The evaluation of our al gorithm on a consecutive series of patients with proximal humeral fractures and, additionally, the quantitative validation on cadaver fractures demonstrated, that the shape of the fracture surface en codes sufficient information to perform the reconstruction with out needing a reconstruction template such as the contralateral anatomy. Further research will focus on the development of strate gies to reduce the number of potential candidates in the pairwise matching at an earlier stage. In addition, it will be interesting to evaluate the application of the method for the anatomical recon struction of fractures of further anatomies, i.e. fractures of the dis tal radius or of the proximal femur.

References

Albrecht, T. , Lüthi, M. , Gerig, T. , Vetter, T. , 2013. Posterior shape models. Med. Image Anal. 17, 959–973 . Albrecht, T. , Vetter, T. , 2012. Automatic fracture reduction. In: Mesh Processing in Medical Image Analysis. Springer, pp. 22–29 . Anakwenze, O.A. , Zoller, S. , Ahmad, C.S. , Levine, W.N. , 2014. Reverse shoulder arthro plasty for acute proximal humerus fractures: a systematic review. J. Shoulder Elbow Surg. 23, e73–e80 . Bicknell, R.T. , DeLude, J.A . , Kedgley, A .E. , Ferreira, L.M. , Dunning, C.E. , King, G.J. , Faber, K.J. , Johnson, J.A. , Drosdowech, D.S. , 2007. Early experience with comput erassisted shoulder hemiarthroplasty for fractures of the proximal humerus: development of a novel technique and an in vitro comparison with traditional methods. J. Shoulder Elbow Surg 16, S117–S125 . Boileau, P. , Krishnan, S.G. , Tinsi, L. , Walch, G. , Coste, J.S. , Mole, D. , 2002. Tuberosity malposition and migration: reasons for poor outcomes after hemiarthroplasty for displaced fractures of the proximal humerus. J. Shoulder Elbow Surg 11, 401–412 . Chowdhury, A.S. , Bhandarkar, S.M. , Robinson, R.W. , Yu, J.C. , 2009. Virtual multi fracture craniofacial reconstruction using computer vision and graph matching. Comput. Med. Imaging Graphics 33, 333–342 . Cuff, D.J. , Pupello, D.R. , 2013. Comparison of hemiarthroplasty and reverse shoulder arthroplasty for the treatment of proximal humeral fractures in elderly patients. J. Bone Joint Surg. Am. 95, 2050–2055 . Cui, M. , Wonka, P. , Razdan, A. , Hu, J. , 2007. A new image registration scheme based on curvature scale space curve matching. Visual Comput. 23, 607–618 . Cvetanovich, G.L. , Chalmers, P.N. , Verma, N.N. , Nicholson, G.P. , Romeo, A .A . , 2016. Open reduction internal fixation has fewer shortterm complications than shoulder arthroplasty for proximal humeral fractures. J. Shoulder Elbow Surg. 25, 624–631 e623 . De Boor, C. , 1978. A Practical Guide to Splines. SpringerVerlag, New York . DeLude, J.A. , Bicknell, R.T. , MacKenzie, G.A. , Ferreira, L.M. , Dunning, C.E. , King, G.J. , Johnson, J.A. , Drosdowech, D.S. , 2007. An anthropometric study of the bilateral anatomy of the humerus. J. Shoulder Elbow Surg. 16, 477–483 . Descoteaux, M. , Audette, M. , Chinzei, K. , Siddiqi, K. , 2006. Bone enhancement filter ing: application to sinus bone segmentation and simulation of pituitary surgery. Comput. Aided Surg. 11, 247–255 . Dijkstra, E.W. , 1959. A note on two problems in connexion with graphs. Numerische Mathematik 1, 269–271 . Fornaro, J. , Keel, M. , Harders, M. , Marincek, B. , Székely, G. , Frauenfelder, T. , 2010. An interactive surgical planning tool for acetabular fractures: initial results. J. Orthopaedic Surg. Res. 5, 50 . Fucentese, S.F. , Sutter, R. , Wolfensperger, F. , Jost, B. , Gerber, C. , 2014. Large metaphy seal volume hemiprostheses for complex fractures of the proximal humerus. J. Shoulder Elbow Surg. 23, 427–433 .

Fürnstahl, P., 2010. Computerassisted planning for orthopedic surgery. Diss., Eid genössische Technische Hochschule ETH Zürich, Nr. 19102, 2010.

Fürnstahl, P. , Szekely, G. , Gerber, C. , Hodler, J. , Snedeker, J.G. , Harders, M. , 2012. Com puter assisted reconstruction of complex proximal humerus fractures for preop erative planning. Med. Image Anal. 16, 704–720 . Gallinet, D. , Clappaz, P. , Garbuio, P. , Tropet, Y. , Obert, L. , 2009. Three or four parts complex proximal humerus fractures: hemiarthroplasty versus reverse prosthe sis: A comparative study of 40 cases. Orthopaedics Traumatol. 95, 48–55 .

L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156 155

Gass, T. , Szekely, G. , Goksel, O. , 2014. Simultaneous segmentation and multiresolu tion nonrigid atlas registration. IEEE Trans. Image Process. 23, 2931–2943 . Gerber, C. , Hersche, O. , Berberat, C. , 1998. The clinical relevance of posttraumatic avascular necrosis of the humeral head. J. Shoulder Elbow Surg. 7, 586–590 . Gerber, C. , Werner, C.M. , Vienne, P. , 2004. Internal fixation of complex fractures of the proximal humerus. J. Bone Joint Surg. Br. 86, 848–855 . Grubhofer, F. , Wieser, K. , Meyer, D.C. , Catanzaro, S. , Beeler, S. , Riede, U. , Gerber, C. , 2016. Reverse total shoulder arthroplasty for acute headsplitting, 3 and 4part fractures of the proximal humerus in the elderly. J. Shoulder Elbow Surg. 25, 1690–1698 . Harders, M. , Barlit, A. , Gerber, C. , Hodler, J. , Székely, G. , 2007. An optimized surgical planning environment for complex proximal humerus fractures. MICCAI Work shop on Interaction in Medical Image Analysis and Visualization . Horn, B.K. , 1987. Closedform solution of absolute orientation using unit quater nions. JOSA A 4, 629–642 . Huffman, G.R. , Itamura, J.M. , McGarry, M.H. , Duong, L. , Gililland, J. , Tibone, J.E. , Lee, T.Q. , 2008. Neer Award 2006: biomechanical assessment of inferior tuberosity placement during hemiarthroplasty for fourpart proximal humeral fractures. J. Shoulder Elbow Surg. 17, 189–196 . Iannotti, J. , Baker, J. , Rodriguez, E. , Brems, J. , Ricchetti, E. , Mesiha, M. , Bryan, J. , 2014. Threedimensional preoperative planning software and a novel informa tion transfer technology improve glenoid component positioning. J. Bone Joint Surg. Am. 96, e71 . JimenezDelgado, J.J. , PaulanoGodino, F. , PulidoRamRamirez, R. , JimenezPerez, J.R. , 2016. Computer assisted preoperative planning of bone fracture reduction: Sim ulation techniques and new trends. Med. Image Anal. 30, 30–45 . Jolliffe, I. , 2002. Principal Component Analysis. Wiley Online Library . Kishon, E. , Wolfson, H. , 1987. 3D curve matching. In: Proceeding of the AAAI Work shop on Spatial Reasoning and Multisensor Fusion, pp. 250–261 . Kovler, I. , Joskowicz, L. , Weil, Y.A. , Khoury, A. , Kronman, A. , Mosheiff, R. , Lieber gall, M. , Salavarrieta, J. , 2015. Haptic computerassisted patientspecific preop erative planning for orthopedic fractures surgery. Int. J. Comput. Assist. Radiol. Surg. 10, 1535–1546 . Kreyszig, E. , 1959. Differential Geometry. The University of Totonto Press, Toronto . Kronman, A. , Joskowicz, L. , 2013. Automatic bone fracture reduction by fracture con tact surface identification and registration, Biomedical Imaging (ISBI). In: 2013 IEEE 10th International Symposium on, pp. 246–249 . Lecouvet, F.E. , Simoni, P. , Koutaïssoff, S. , Vande Berg, B.C. , Malghem, J. , Dubuc, J.E. , 2008. Multidetector spiral CT arthrography of the shoulder: Clinical applications and limits, with MR arthrography and arthroscopic correlations. Eur. J. Radiol. 68, 120–136 . Levy, J.C. , Everding, N.G. , Frankle, M.A. , Keppler, L.J. , 2014. Accuracy of patientspe cific guided glenoid baseplate positioning for reverse shoulder arthroplasty. J. Shoulder Elbow Surg. 23, 1563–1567 . Lorensen, W.E. , Cline, H.E. , 1987. Marching cubes: a high resolution 3D surface con struction algorithm. In: ACM Siggraph Computer Graphics. ACM, pp. 163–169 .

Lüthi, M., Jud, C., Gerig, T., Vetter, T., 2016. Gaussian Process Morphable Models. arXiv preprint arXiv:1603.07254 . MATLAB, 2015. version 8.5 (R2015a). The MathWorks (Inc.). McBride, J.C. , Kimia, B.B. , 2003. Archaeological fragment reconstruction using curve– matching. In: 2003 Conference on Computer Vision and Pattern Recognition Workshop, p. 3 . Moghari, M.H. , Abolmaesumi, P. , 2008. Global registration of multiple bone frag ments using statistical atlas models: feasibility experiments. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2008, 5374–5377 .

Murase, T. , Oka, K. , Moritomo, H. , Goto, A. , Yoshikawa, H. , Sugamoto, K. , 2008. Threedimensional corrective osteotomy of malunited fractures of the upper ex tremity with use of a computer simulation system. J. Bone Joint Surg. Am. 90, 2375–2389 . Neer, C.S.,2nd , 1970. Displaced proximal humeral fractures. II. Treatment of three– part and fourpart displacement. J. Bone Joint Surg. Am. 52, 1090–1103 . Nguyen, D. , Ferreira, L.M. , Brownhill, J.R. , King, G.J. , Drosdowech, D.S. , Faber, K.J. , Johnson, J.A. , 2009. Improved accuracy of computer assisted glenoid implanta tion in total shoulder arthroplasty: an invitro randomized controlled trial. J. Shoulder Elbow Surg. 18, 907–914 . Okada, T. , Iwasaki, Y. , Koyama, T. , Sugano, N. , Chen, Y.W. , Yonenobu, K. , Sato, Y. , 2009. Computerassisted preoperative planning for reduction of proximal femoral fracture using 3DCT data. IEEE Trans. Biomed. Eng. 56, 749–759 . Olsson, P. , Nysjö, F. , Hirsch, J.M. , Carlbom, I.B. , 2013. A hapticsassisted craniomax illofacial surgery planning system for restoring skeletal anatomy in complex trauma cases. Int. J. Comput. Assist. Radiol. Surg. 8, 887–894 . Papaioannou, G. , Theoharis, T. , 2003. Fast fragment assemblage using boundary line and surface matching. Computer Vision and Pattern Recognition Workshop, 2003. CVPRW ’03. Conference on pp. 2–2 . Salomon, D. , 2007. Curves and Surfaces for Computer Graphics. Springer Science & Business Media . Shukla, D.R. , McAnany, S. , Kim, J. , Overley, S. , Parsons, B.O. , 2016. Hemiarthroplasty versus reverse shoulder arthroplasty for treatment of proximal humeral frac tures: a metaanalysis. J. Shoulder Elbow Surg. 25, 330–340 . Skedros, J.G. , Knight, A.N. , Pitts, T.C. , O’Rourke, P.J. , Burkhead, W.Z. , 2016. Radio graphic morphometry and densitometry predict strength of cadaveric proxi mal humeri more reliably than age and DXA scan density. J. Orthop. Res. 34, 331–341 . Thomas, T.P. , Anderson, D.D. , Willis, A.R. , Liu, P. , Frank, M.C. , Marsh, J.L. , Brown, T.D. , 2011. A computational/experimental platform for investigating threedimen sional puzzle solving of comminuted articular fractures. Comput. Methods Biomech. Biomed. Eng. 14, 263–270 . Üçoluk, G. , Toroslu, I.H. , 1999. Automatic reconstruction of broken 3D surface ob jects. Comput. Graphics 23, 573–582 . Umeyama, S. , 1991. Leastsquares estimation of transformation parameters between two point patterns. IEEE Trans. Pattern Anal. Mach. Intell. 376–380 . Vlachopoulos, L. , Dünner, C. , Gass, T. , Graf, M. , Goksel, O. , Gerber, C. , Székely, G. , Fürnstahl, P. , 2016a. Computer algorithms for threedimensional measurement of humeral anatomy: analysis of 140 paired humeri. J. Shoulder Elbow Surg. 25 e38–e48 . Vlachopoulos, L. , Schweizer, A. , Meyer, D.C. , Gerber, C. , Furnstahl, P. , 2016b. Threedi mensional corrective osteotomies of complex malunited humeral fractures using patientspecific guides. J. Shoulder Elbow Surg. 25, 2040–2047 . Wijgman, A.J. , Roolker, W. , Patt, T.W. , Raaymakers, E.L. , Marti, R.K. , 2002. Open re duction and internal fixation of three and fourpart fractures of the proximal part of the humerus. J Bone Joint Surg Am 84a, 1919–1925 .

Willis, A., Anderson, D., Thomas, T., Brown, T., Marsh, J.L., 2007. 3D reconstruction of highly fragmented bone fractures, pp. 65121P65121P65110.

Winkelbach, S. , Rilk, M. , Schönfelder, C. , Wahl, F.M. , 2004. Fast Random Sample Matching of 3d Fragments. In: Rasmussen, C.E., Bülthoff, H.H., Schölkopf, B., Giese, M.A. (Eds.), Pattern Recognition: 26th DAGM Symposium, Tübingen, Ger many, August 30 September 1, 2004. Proceedings. Berlin Heidelberg, Berlin, Heidelberg. Springer, pp. 129–136 . Winkelbach, S. , Wahl, F.M. , 2008. Pairwise matching of 3D fragments using cluster trees. Int. J. Comput. Vision 78, 1–13 . Zhou, B. , Willis, A. , Sui, Y. , Anderson, D.D. , Brown, T.D. , Thomas, T.P. , 2009. Virtual 3D bone fracture reconstruction via interfragmentary surface alignment. In: 2009 IEEE 12th International Conference on Computer Vision Workshops, ICCV Work shops 2009, pp. 1809–1816 .

156 L. Vlachopoulos et al. / Medical Image Analysis 43 (2018) 142–156

Lazaros Vlachopoulos received his medical degree in 2004 from the RWTH Aachen, Germany, and passed in 2011 the Swiss board exams in orthopaedic surgery. In 2017 he received his PhD in medical image analysis from the ETH Zurich, Switzerland. He is currently working as a consultant orthopaedic surgeon at the Balgrist University Hospital in Zurich, Switzerland, specialized in computerassisted orthopaedic surgery.

Gábor Székely received the Graduate degree in chemical engineering, the Graduate degree in applied mathematics, and the PhD degree in analytical chemistry from the Technical University of Budapest and Eötvös Lórand University, Budapest, Hungary, in 1974, 1981, and 1985, respectively. He has been working at the Computer Vision Laboratory, ETH Zurich, Switzerland, as a Full Professor for Medical Image Analysis and Visualization, where he has been involved in the development of image analysis, visualization, and simulation methods for computer support of medical diagnosis, therapy, training, and education.

Christian Gerber graduated 1977 from the medical school at the University of Berne, Switzerland. He was trained in orthopaedic surgery and specialized in shoulder surgery at the University of Texas, San Antonio in 1984 and subsequently trained on tumor foot and ankle surgery as well as pediatric orthopaedics in Paris in 1985. In 1991 he was promoted to associate professor in orthopaedic surgery. In 1992, he was appointed chief of the Department of Orthopaedics in Fribourg, Switzerland, and since 1995 he is at the Balgrist University Hospital in Zurich, Switzerland, where he holds the position of medical director and Chairman of the Department of Orthopaedics.

Philipp Fürnstahl received the MSc degree in technical mathematics and information procession from the Technical University of Graz, Austria, in 2005. In 2010, he received the PhD degree in medical image analysis from the ETH Zurich, Switzerland, for his research in the field of computerassisted preoperative surgery planning. Philipp Fürnstahl is currently head of the Computer Assisted Research and Development (CARD) Group at the University Hospital Balgrist in Zurich, Switzerland, where he is involved in the development of patientspecific solutions for orthopaedic surgeries.