Int J CARS DOI 10.1007/s11548-014-1011-2

ORIGINAL ARTICLE

Computer-assisted fracture reduction: a new approach for repositioning femoral fractures and planning reduction paths Jan Buschbaum · Rainer Fremd · Tim Pohlemann · Alexander Kristen

Received: 10 January 2014 / Accepted: 16 April 2014 © CARS 2014

Abstract Purpose Reduction is a crucial step in the surgical treatment of bone fractures to achieve anatomical alignment and facilitate healing. Surgical planning for treatment of simple femoral fractures requires suitable gentle reduction paths. A plan with optimal movement of fracture fragments from the initial to the desired target position should improve the reduction procedure. A virtual environment which repositions the fracture fragments automatically and provides the ability to plan reduction paths was developed and tested. Methods Virtual 3D osseous fragments are created from CT scans. Based on the computed surface curvatures, strongly curved edges are selected and fracture lines are generated. After assignment of matching points, the lines are compared and the desired target position is calculated. Planning of reduction paths was achieved using a reference-coordinatesystem for the computation of reduction parameters. The fracture is reduced by changing the reduction parameters step by step until the target position is reached. To test this system, nine different fractured SYNBONE models and one human fracture were reduced, based on CT scans with varying resolution. Results The highest mean translational error is 1.2 ± 0.9 (mm), and the rotational error is 2.6 ± 2.8 (◦ ), both of which are considered as clinically acceptable. The reduction paths J. Buschbaum (B) · R. Fremd Fachbereich Angewandte Ingenieurwissenschaften, Fachhochschule Kaiserslautern - University of Applied Sciences, Morlauterer Straße 31, 67657 Kaiserslautern, Germany e-mail: [email protected] A. Kristen · T. Pohlemann Klinik für Unfall-, Hand- und Wiederherstellungschirurgie der Universitätsklinik des Saarlandes, Homburg/Saar, Geb. 57 Kirrberger Straße 1, 66421 Homburg, Saar, Germany

can be planned manually or semi-automatically for each fracture. Conclusions Automated fracture reduction was achieved using a system based on preoperative CT scans. The automated system provides a clinically feasible basis for planning optimal reduction paths that may be augmented by further computer- or robot-assisted applications. Keywords Computer-assisted fracture reduction · Femoral fractures · Virtual 3D models alignment · Fitting fracture lines · Reduction parameters · Planning reduction paths Introduction The treatment of femoral fractures is a common task in orthopedic surgery. The procedure is performed in two main steps. First, the fracture is reduced, and then, the bone fragments are fixed by plates or intramedullary nails. A crucial step is the reduction, whose goal is to reposition the bones in their correct anatomical alignment [1]. A common problem is to deduce the desired target position from medical imaging. This can lead to postoperative malalignment which has great impact on the course of healing. For increasing the accuracy and reducing the risk of postoperative malalignment, a variety of computer-aided navigation systems were developed [2–11]. Further difficulties are high forces which occur during fracture reduction. These high forces increase the physical load of the surgeons and can prevent reduction movements [12,13]. To overcome these forces, the developments range from simple auxiliary devices [14–16] through to robotassisted systems [17–29]. In this field, only a few studies have investigated in planning reduction paths for an automatic reduction procedure

123

Int J CARS

[13,20,30]. In particular, we see these pathways as the key to success of a facile and gentle reduction procedure [31]. The goal of our research activities is to include muscle strength, soft tissue deformations and bone collisions for planning of optimal reduction paths. These paths are characterized by an optimal sequence of motion patterns to accurately and carefully reduce the fragments. Unnecessary and prohibited movements that can lead to high forces, damage muscle, tendons or other vulnerable structures are to be avoided. The so planned paths can be used for navigated fracture reduction. In this case, the surgeon will be guided by a navigation device along the paths from the initial to the target position. For this purpose, we developed a virtual environment which automatically repositions the fracture and especially allows the planning of reduction paths. This is described in the method of this paper. It should be noted that we only work with simple femoral shaft fractures. An extension to complex fractures or other fracture types is the task of further work. Methods The presented method is divided into two parts. The first part is the automatic reduction algorithm. The bone fragments are segmented from a CT scan, and virtual 3D models are created. Based on the computed surface curvatures, strongly curved break lines can be distinguished from the smooth surface of the bone. In the areas of high convex curvature crest points are extracted and connected to crest lines, denoted as fracture lines. By using special surface properties, the fracture edges are reconstructed and the fragments are repositioned. In the second part, we present a path-planning approach. By utilizing a reference-coordinate-system, the computation of special reduction parameters is possible. These describe the displacement and malalignment in the current position and are used for planning manual and semiautomatic reduction paths. For this, both fragments are visualized on a screen and the fracture is reduced by changing the reduction parameters step by step until the desired target position is reached. Automatic reduction in bone fractures through reconstruction of fracture lines Serving as a basis for the method, the distal (Dist) and the proximal (Pr ox) fragment are segmented out of a CT scan. By extracting isosurfaces with the same gray values (Hounsfield units), virtual 3D models are created [32]. These 3D models are represented by the point set P := {P1 , P2 , . . . , Pn } and the triangular surfaces F := {Fi |i = 1, 2, . . . , m}. Additionally, the surface normal vectors   S := S j | j = 1, 2, . . . , n for each point are computed [33].

123

The goal of the reduction method is to compute a rigid transformation T which aligns the distal fragment to the proximal fragment. Therefore, fracture lines are generated on triangular meshes. By using surface curvature properties, an estimated transformation is assumed. After an assessment and assigning of matching fracture lines, the desired transformation is calculated and the fragments are repositioned. Extraction of fracture lines on triangular meshes Fracture lines are special crest lines (or ridge lines) in regions with high convex surface curvature. To extract fracture lines on triangular meshes, the computation of the surface curvature properties is needed. After aligning the point set to the surface normal vectors, a function φ(u, v) = au 2 + bv 2 + cuv +du +ev + f is fitted in the neighborhood of point P by the method of least squares [34]. The Weingarten endomorphism L (or shape operator) is built from the partial derivatives of that function φ: ⎞ ⎛ 2 2 ∂ φ 2u

∂ L=⎝

∂ φ ∂u∂v

∂2φ ∂2φ ∂v∂u ∂ 2 v

⎠

(1)

The eigenvalues of L (Eq. 1) are the principal curvatures k1 and k2 . It is |k1 | > |k2 |. The corresponding eigenvectors of L are the principal curvature directions t1 and t2 [35,36]. Therefore, t1 is in the direction of greatest curvature k1 . In summary, for each point P j ∈ P of the triangular mesh, the surface curvatures k1 (P j ), k2 (P j ) and the curvature directions t1 (P j ), t2 (P j ) are computed. Figure 1a, d shows a heat map of the computed surface curvature. The fracture edges (drawn in blue) can be distinguished from other areas by mean curvature H in Eq. (2). For the strongly convex curved edges applies [33]: H (P j ) =

k1 (P j ) + k2 (P j ) < 0 ⇒ k1 (P j ) < 0, 2 because |k1 (P j )| > |k2 (P j )|

(2)

The selected points of the fracture edge are B := {P ∈ P|k1 (P) < ε < 0}. Thereby ε is a threshold to specify the magnitude of the convex surface curvature. In order to generate the fracture lines, the crest points must be extracted out of B. As described in [37,38], crest points are local extrema of curvature. They are defined as the zero crossing of the curvature derivative in their corresponding directions: ei = ∇ki · ti = 0, i = 1, 2

(3)

The values ei differ due to the curvature direction and due to the convex or concave curved surface. The required crest points are maxima of convex surfaces in the first principal curvature direction. Therefore, Eq. (3) becomes e1 = ∇k1 · t1 = 0.

Int J CARS Fig. 1 a Heat map of surface curvature of the proximal fragment (P r ox). The high convex curvatures of the fracture edge are drawn in blue. b Extracted crest points of P r ox with the second curvature direction t2 . c Generated fracture lines F L P of P r ox. d Heat map of surfaces curvature of the distal fragment (Dist). The high convex curvatures of the fracture edge are drawn in blue. e Extracted crest points of Dist with the second curvature direction t2 . f Generated fracture lines F L D of Dist

The calculation of e1 on discrete triangular mesh is performed as described in [39]. We denoted P j+1 as the nearest point of the neighborhood of P j ∈ B which lies in the positive t1 direction. P j−1 is the same only in negative t1 direction. The curvatures difference (Eq. 4) between the points can be regarded as an approximation of e1 [39]: dk j+1 = k1 (P j+1 ) − k1 (P j ) and dk j−1 = k1 (P j−1 ) −k1 (P j )

(4)

Taking the sign (k1 < 0) into account, point P j is a convex crest point if dk j+1 > 0 and dk j−1 > 0. The crest points of the fracture edge are CP := {P ∈ B|(dk j+1 > 0) ∧ (dk j−1 > 0)} and are shown in Fig. 1b, e. Due to the surfaces properties, the generation of the fracture lines out of CP is possible in a simple manner. As seen in Fig. 1b, e, the surface curvature direction t2 of CP specifies the curve shapes of the fracture lines. By tracing the CP along their t2 direction, crest points can be connected as crest

lines [40]. Therefore, the nearest crest point C Pi+1 ∈ CP in the positive t2 direction of the neighborhood is selected and connected with C Pi ∈ CP by a line. For C Pi−1 ∈ CP, the same applies for the negative t2 direction. In a dense order of points, this leads to jagged lines. Additional criteria, such as a minimum distance d and a permissible angular deviation ε improved the results. By d, the length of the each line segment can be adjusted and ε controls the roughness of the line. Out of a large amount of crest points, these special fracture lines of the distal FL D (Fig. 1c) and the proximal fragment FL P (Fig. 1f) are created and serve as the basis of the following method. Computation of the target position by matching fracture lines The goal is to compute a rigid transformation, whereby the opposite fracture lines FL D and FL P are mapped together (Eq. 5).

123

Int J CARS

Distal

Proximal

Fig. 2 Geometrical relationship of fracture lines. The points of the fracture lines are drawn in blue. Additionally, the surface frame ( S pi , t1 pi , t2 pi ) at point pi on the proximal fragment and a corresponding frame ( Sd j , t1d j , t2d j ) on the opposite distal fragment at point di are shown

T : FL D → FL P ⇔ T (FL D ) = FL P

(5)

For registration of those point clouds, the iterative closest point algorithm (ICP) [41] is a commonly used method. The spatially nearest points are assigned to each other, and a transformation is calculated. The algorithm iterates until the sum of the squared distances reaches a minimum. The problem is that the nearest points are not necessarily matching points. Incorrect assignments of the corresponding points and possible outliers have a fatal influence on the result of the ICP. The main problem particularly lies in finding these matching points. Invariant parameters such as curvature and torsion are often used for comparing curves and finding geometrically similar lines [42]. But the numerical calculation leads to noisy values. Only through an additional smoothing method [43] or the generation of special parameters [44,45], stable curvature features can be calculated. However, line segments from the inner and outer tube of the bone can have similar characteristics whereby an unambiguous assignment is difficult. As described in [42], the surface properties can be used for a registration method. The surface frames consist of the surface normal vector, first and second surface curvature direction vector and can be regarded as the Frenet frame of a curve. For our approach, we can use surface frames by taking the geometric relationships into account. These are shown in Fig. 2. Additionally, to the point location pi ∈ FL P on the proximal fragment and dj ∈ FL D on the distal site, their corresponding surface frames F pi and Fd j are drawn. In Eq. (6), these are represented in homogeneous coordinates. This notation is commonly used in robotic applications [46].  F pi :=



t2 pi S pi t1 p j pi ; 0 0 0 1

123

 Fd j :=



t2d j Sd j t1d j dj ; 0 0 0 1 (6)

It can be seen that Fd j is rotated by ϕ = π2 around the t2 -axis with respect to F pi . Taking this geometric relationship into account, the transformation which aligns the fragments can be calculated as shown in Eq. (7). ϕ

· Fd j Td j pi = F p−1 i

(7)

ϕ

For that Fd j is the frame rotated by ϕ =

π 2

around the t2 -axis.

The new alignment of distal points d ∈ FL D , after per How well the forming the transformation, is d∗ = Tdi pi · d. fracture lines fit together is assessed by the number of best matching points (bmp) in Eq. (8): 

(8) bmp d∗ , FL P := arg min d∗ − pi < εdist ∀ pi ∈F L P

For all points pi ∈ FL P , the point with the smallest Euclidean distance in d∗ ∈ FL D is searched for and assigned. If the distance is smaller than the threshold εdist , both point pairs are stored as bmp. The best matching transformation is when the maximum number of fracture lines is fitted that is if many bmp are found. The aim is to find that point pair pi and dj and thus Td j pi whereby most of the fracture lines are matched. This transformation is only a fair approximation of the desired target position. The exact transformation can be calculated from the matching points. We denote X ∈ FL P and Y ∈ FL D as the matching points bmp of the fracture lines from Eq. (8). The two point sets xi ∈ X and yi ∈ Y are displaced by t and rotated by R so that the sum of the squared distances in Eq. (9) is minimized: E=

n

(xi − R · yi − t)2

(9)

i=1

To solve R and t, we use the singular value decomposition (svd) method as described in [47].

Int J CARS

Fig. 3 Automatic reduced fracture of a SYNBONE model (fracture type 32-A2 classified according to AO rules [1]). Three different views and the corresponding points (bmp) are shown

As the result, we compute a 4 × 4 transformation matrix T (Eq. 10) which fits the fracture lines together and solves Eq. (5).  t3×1 R3×3 (10) T := 0 0 0 1 In principle, all points pi ∈ FL P and d∗ ∈ FL D can be tested, but there can be a variety of matches that lead to success. Our method uses a random-based approach [48] which greatly reduces the computational effort. Therefore, pi ∈ FL P and dj ∈ FL D are randomly selected points and Tdi pi is the calculated transformation. The algorithm terminates when an optimum number of best matching points bmp > εmatch is reached. Otherwise, another randomly chosen point pair is tested until bmp > εmatch is achieved or all points have been checked. Finally, the distal fragment is repositioned by this transformation from the initial position Distinit to the desired target position: Disttarget = T · Distinit

(11)

The result of a reduced fracture is shown in Fig. 3. Approach for planning reduction paths We define the reduction path as the trajectory ΓDist of the distal fragment from the initial to the target position: Disttarget = ΓDist · Distinit

(12)

The path is a composition of the single movements which are required for assembling the fracture. With respect to Eq. (13), T˜i are the transformations which describe the alignment in each reduction step i. The combination of all T˜i leads to the reduction path ΓDist by which the distal fragment is moved from the initial (i = 1) to the target position (i = N ): ΓDist :=

1

T˜i = T˜N · T˜N −1 . . . · T˜1

(13)

i=N

To plan such a trajectory, we utilized a reference-coordinatesystem and calculated special reduction parameters. Then the fractures can be reduced manually or semi-automatically. Calculation of reduction parameters The procedure to generate a reference-coordinate-system is shown in Fig. 4. After the repositioning of the distal fragment (Fig. 4a), a target coordinate system Ctarget is generated in the region of the fracture (Fig. 4b). The main axis of the coordinate system is an approximation of the middle axis of the bone shaft. For this, it is essential that a high number of surrounding vertex normal vectors are orthogonally oriented. Like described in [48], we create various axis-hypotheses H (Eq. 14) by the surface normal vectors of two randomly selected points Si and S j : H :=

Si × S j Si × S j

(14)

123

Int J CARS

Fig. 4 Procedure to generate target and initial coordinate systems. a The initial position of the fragments. b The repositioned fracture and the computed target coordinate system Ctarget . c Visualization of Ctarget and by T −1 retransformed initial coordinate system Cinit Table 1 Description and calculation of the six reduction parameters. T is the relative transformation from Eq. (17)

Parameter Description

Calculation

ds (m)

Ventral/dorsal displacement (sagittal direction)

ds := T1,4

dt (m)

Lateral/medial displacement (transversal direction)

dt := T2,4

dl (m)

Distal/proximal displacement (longitudinal direction) dl := T3,4

 T T1,1 Inner/outer rotational malalignment α = atan2 cos2,1 β , cos β  

2 + T 2 Antecurvation/retrocurvation malalignment β = atan2 −T3,1 , T1,1 2,1 

T T3,3 Varus/valgus malalignment γ = atan2 cos3,2 β , cos β

α (◦ ) β (◦ ) γ (◦ )

The hypotheses are evaluated by checking how many of the total normal vectors are in orthogonal direction to the axis. The best rating hypotheses is selected as the fracture middle axis and denoted as l (longitudinal direction). The other two axes of the coordinate system are orthogonal to l and oriented in transversal t and sagittal s directions. The origin o of the coordinate system is the centroid of the fracture lines xi ∈ X . The complete Cartesian target coordinate system is: 

Ctarget

s t l o := 0 0 0 1

(15)

By the inverse matrix of T , the target coordinate system is transformed back into the initial position: Cinit = T −1 · Ctarget

(16)

The so generated coordinate systems Cinit and Ctarget are shown in Fig. 4c. With the relative transformation matrix in Eq. (17) between both coordinate systems, six reduction parameters are calculated [46] and listed in Table 1. −1 · Cinit T = Ctarget

(17)

These reduction parameters are used for the description of the fragment position and provide information about the exact translational and rotational deviations of the distal fragment from the target position. Thus, the basis for planning paths is created.

123

Manual and semi-automatic planning of reduction paths For manual or semi-automatic creation of paths, both fragments are visualized on a screen. By changing the reduction parameters, the fracture is reduced. Manual reduction means that by changing each single parameter, the distal fragment is moved step by step in the desired target position. The semi-automatic approach uses special movement patterns, which are listed in Table 2. The illustrated example in Fig. 5 shows a sequence of five simple reduction movements by which the fracture is reduced. Starting from the initial position, the middle axes are aligned parallel before a distraction in distal direction is performed. Subsequently, the both middle axes are fitted together. The following rotation around the middle axis and the merging completes the reduction process. The whole reduction maneuver is described by the reduction parameters and is listed in Table 2. The order of movements can be varied, and thus, different paths can be created. The sequence of motion is connected to a reduction path related to Eq. (13). For each reduction step (1 . . . N ), the relative transformation is built: ⎡ ⎤ ds ⎢ R(α,β,γ ) dt ⎥ ⎥ (18) T1...N := ⎢ ⎣ dt ⎦ 0 0 0 1 Whereby R(α,β,γ ) = Rl(α) · Rt(β) · Rs (γ ) is the rotation matrix with the ZYX Euler angles around the coordinate axes. The current pose of the fracture-coordinate-system is: C1...N = Ctarget · T1...N

(19)

Int J CARS Table 2 List of the individual reduction steps as seen in Fig. 5 with their corresponding reduction parameters Step

Movement pattern

Reduction parameter ds (m)

dt (m)

dl (m)

α (◦ )

β (◦ )

γ (◦ )

0

Initial position

0

0.04

−0.055

14

−1.75

6

1

Align middle axis

0

0.04

−0.055

14

0

0

2

Distraction

0

0.04

0.01

14

0

0

3

Join middle axis together

0

0

0.01

14

0

0

4

Adjust rotational malalginment

0

0

0.01

0

0

0

5

Join together (target position)

0

0

0

0

0

0

Fig. 5 A sequence of individual reduction movements conducted as a reduction path. a Step 0: the initial pose of both fragments. b Step 1: align the middle axis. c Step 2: distraction in longitudinal direction.

d Step 3: join middle axis together. e Step 4: adjust inner or outer rotational malalignment. f Step 5: join bones together (target position)

With C1...N from Eq. (19) the single transformations for each step T˜1...N are:

The result of a reduced SYNBONE fracture (case A, type 32-A2) is shown in Fig. 3. Table 3 shows the MTEs and SDs of six selected and different fractures. The left column classifies the CT Resolution (A, B, C) and the different fracture types according to the AO rules [1]. The six rightmost columns indicate the MTEs and SDs for the ventral/dorsal, lateral/medial, distal/proximal displacement and inner/outer, antecurvation/retrocurvation, varus/valgus malalignment. The highest values for translational deviations are 1.2 ± 0.9 (mm) and for rotational errors 2.6 ± 2.8 (◦ ), which are considered as clinically acceptable. In addition, the results show that the method can be applied to various fracture types of human cadaveric bones or SYNBONES models with different CT resolutions. It should be noted, if the fragments are segmented correctly, it does not matter whether it is human bone or SYNBONE model. As a further result, we present a method which allows planning reduction paths for each fracture. As a new element, reduction parameters are introduced. These describe the translational and rotational deviations of the distal fragment from its target position. By adjusting this displacement and malalignment, the two fragments are precisely fitted together. This can be done manually or semi-automatically. By varying the sequences of movements or movement patterns, paths can easily be planned. Thus, it is possible to plan different paths for each fracture. Such a created path is shown

T˜1 T˜2 .. . T˜N

−1 = C1 · Cinit = C2 · C1−1 .. .. . . = C N · C N−1−1

(20)

Finally, we can connect the single transformation to a reduction path ΓDist as shown in Eq. (13).

Results For testing our method, we fractured 9 SYNBONE models (type 2162 SYNBONE AG, Malans, Switzerland) in different ways and recorded several CT scans with different resolutions. The voxel size of case A is 0.58 mm × 0.58 mm × 0.5 mm and case B is 0.85 mm×0.85 mm×2.0 mm. In addition, the method was tested on a human cadaveric fracture (case C with a resolution of 0.34 mm × 0.34 mm × 1.0 mm). The evaluation of the reduction algorithm is difficult because of the missing ground truth. For that, we used a statistical test. Each fracture was reduced 100 times. The results were checked visually, and the mean target errors (MTE) and the standard deviations (SDs) were measured.

123

Int J CARS Table 3 Mean target errors (MTE) and standard deviations (SD) of the automatic reduction algorithm

CT Resolution and fracture type

MTE ± SD ds (mm)

dt (mm)

dl (mm)

α (◦ )

β (◦ )

γ (◦ )

32-A2

0.1 ± 0.1

0.9 ± 0.5

0.2 ± 0.1

2.1 ± 1.5

0.7 ± 0.4

1.9 ± 1.5

32-A3

0.2 ± 0.1

0.1 ± 0.1

0.2 ± 0.2

0.6 ± 0.6

1.8 ± 1.4

0.8 ± 0.7

31-A1

0.2 ± 0.1

0.5 ± 0.4

0.2 ± 0.3

2.5 ± 1.5

2.0 ± 1.8

0.7 ± 0.8

32-A1

0.8 ± 0.6

0.2 ± 0.3

0.7 ± 0.6

1.2 ± 1.0

1.0 ± 0.9

0.8 ± 0.7

33-A1

0.9 ± 0.6

0.2 ± 0.4

1.2 ± 0.9

1.9 ± 2.3

2.6 ± 2.8

1.3 ± 1.1

0.8 ± 0.6

0.9 ± 0.5

0.7 ± 0.4

0.3 ± 0.2

0.1 ± 0.1

1.1 ± 0.9

A

B

C 32-A2

in Table 2 and Fig. 5. It describes an ordered series of movements which are required to move from the initial position to the exact target position. This allows analyzing and planning the tactical procedure of a fracture reduction preoperatively and provides a good basis for developing further planning strategies.

Discussion Reduction is a crucial step in the surgical treatment of bone fractures to achieve anatomical alignment and facilitate healing. A common problem is to deduce the desired target position from medical imaging. This can lead to postoperative malalignment which has great impact on the course of healing [1]. During fracture reduction procedure, further difficulties can arise. High forces which are caused by stretching the surrounding soft tissues increase the physical load of the surgeons and can prevent the desired movement [12,13]. But a further excessive stretching can damage muscles, tendons or other vulnerable structures. An additional obstacle occurs when the bones collide. Then, the movement is obstructed, and it is impossible to achieve the desired position. Finally, when any obstacle occurs, a different or an additional movement has to be performed. This may lead to a trial and error procedure. The consequences can be higher radiation exposure and longer operation times. Therefore, the objective of our research activities is the planning of optimal reduction paths. Through an optimal sequence of motions, the fracture should be accurately and carefully reduced. Unnecessary and prohibited movements, which can lead to collisions and high forces, should be avoided. Particularly, we see these pathways as the key to success of a facile and gentle reduction procedure. The so planned paths can be used for navigated fracture reduction. In this case, the surgeon will be guided by a navigation device along the paths from the initial to the target position.

123

So far, only a few studies have dealt with planning reduction paths. Westphal et al. [20] developed an extensive approach for automated robot-assisted reduction. By probabilistic motion-planning typical reduction movements, denoted as skill primitives, are connected to a path. Based on 3D imaging, the objective of their planning strategy is to minimize the distractions whereby forces should be reduced. This raises the question whether it is sufficient to consider only the distractions for the planning strategy. For simulating forces during fracture reduction, Graham et al. [13] and Joung et al. [30] developed a musculoskeletal model. The results show that the forces can vary greatly due to different trajectories. The influence of the muscle forces on each movement and the outcome for the planning of reduction paths are not described. It will be a task of our further work to plan reduction paths automatically and under the consideration of collisions and simulated forces. For this purpose, we have developed the presented planning approach. The advantage of our method is that the target position is computed first. Through the knowledge of the desired target position, a special referencecoordinate-system can be utilized. This allows the computation of six reduction parameters which describe the displacement and malalignment between the initial and the target position. Thus, a description of reduction movement with respect to the reference-coordinate-system is possible. By manual movements or semi-automatic movement patterns, the fracture can be reduced step by step in a simple way. For this, both fragments are visualized on a screen and the fracture is reduced by changing the reduction parameters until the exact target position is reached. By varying the sequences of movements or movement patterns, different paths can be planned manually or semi-automatically. This allows analyzing and planning the tactical procedure of a fracture reduction and finally the choice of an optimal path. Furthermore, the presented approach provides a good basis for developing further planning strategies.

Int J CARS

One main part of our method is an automatic reduction algorithm which computes the desired target position. A few methods for automatic reduction in virtual 3D femur fractures are known. Winkelbach and Wahl et al. developed several approaches. In [49], the center lines of the cylindrical fragments are detected with a kind of Hough transform. After selecting the fracture surface, they were registered by a Z-buffer method. This approach works on cylindrical fragments, but can fail at fractures near joints. Another approach is a RANSAC method [50], whose aim is to maximize the contact area of the two fragments. A fragment pose is assumed and assessed by the number of points which are in tangential contact. Problems occur when the selected contact area does not correspond to the fracture surface. An improvement is presented in [48]. With assistance of an estimated center line, the points of the fracture surface are selected and only these points are used for the contact assessment. Here too, the required center line and the selection of the fracture surface can lead to problems. By the method of Josokowiscz and Kronman [51], the fracture surfaces are selected based on the curvature of the virtual models and the intensity of the gray values of the CT image. The fracture surfaces are aligned by a principal component analysis and are registered by the ICP algorithm. A disadvantage of ICP is that incorrect assignments of the corresponding surfaces and possible outliers have a fatal influence on the result. Unlike the described methods, the approaches [52,53] reduce the fragments with respect to the intact bone. As a reference either a statistical atlas [52] or the intact contralateral bone [53] are used. Our method reduces the fracture automatically by reconstructing the fracture lines of the fracture surface. By computing the surface curvatures of the triangulated models, which we create from CT scans, crest lines can be generated in the area of highly curved fracture edges. These fracture lines consist of a few particular characteristic points, with geometrical features, which represent almost completely the contour of the fracture edges. By using the properties of the surfaces at the fracture lines, matching points are assigned and subsequently used to calculate the exact transformation. Thereby the two fragments are repositioned. The generation of additional parameters, such as curvature, torsion or special parameters [44,45], are not needed for this method. Due to the use of the fracture lines, this method does not require any symmetry properties, center lines, contra-lateral bone or special a-priori knowledge. Thus, the method is suitable for various femoral fracture types of human bones or SYNBONES models, and it is largely independent of the resolution of the CT recordings. It should be noted, that if the fragments are segmented correctly, it does not matter whether it is human bone or SYNBONE model. So far we have only considered simple femoral shaft fractures, but an extension

to complex fractures or other fracture types should be possible. A further task will be the consideration of additional human fractures to evaluate our method. As a disadvantage of the method, the generation of fracture lines has to be mentioned. This is done by a simple algorithm which can cause that parts of the fracture edge are not detected. By improving this, it should be possible to get better, finer or smoother lines which should have an impact on the overall result.

Conclusion In this paper, we present a virtual environment which automatically reduces a bone fracture and allows the planning of reduction paths. The bone fragments are segmented from a CT scan, and virtual 3D models are created. Special fracture lines are generated based on the computed surface curvatures of the models. By using special surface properties, the fracture edges are reconstructed and the fragments are repositioned. By utilizing a reference-coordinate-system, the computation of special reduction parameters and thus the planning of reduction paths is possible. For this, both fragments are visualized on a screen and the fracture is reduced by changing the reduction parameters until the desired target position is reached. By varying the sequences of movements or movement patterns, different paths can be planned manually or semi-automatically. This allows analyzing and planning the tactical procedure of a fracture reduction preoperatively. The approach provides a good basis for developing strategies for planning optimal reduction paths. This allows the development of further computer- and robot-assisted applications. Additionally, it can deliver fundamental knowledge for the general procedure of fracture reduction. Acknowledgments This research project is a cooperative project between the Department of Trauma, Hand and Reconstructive Surgery of the University Hospital of the Saarland and the University of Applied Sciences Kaiserslautern. Our thanks go to both partners, for technical and financial support. The research project is currently funded by the “Stiftung Rheinland-Pfalz für Innovation.” Conflict of interest Jan Buschbaum, Rainer Fremd, Tim Pohlemann and Alexander Kristen declare that they have no conflict of interest.

References 1. Rüedi T, Buckley RE, Morgan CG (2007) AO principles of fracture management, books and DVD, 2nd edn. Thieme, AO Pub, Stuttgart, New York 2. Gardner MJ, Citak M, Citak M et al (2008) Navigated femoral anteversion measurements: a new intraoperative technique. Injury 39(4):467–471

123

Int J CARS 3. G Zheng, X Dong, X Zhang et al (2005) Automated detection and segmentation of diaphyseal bone fragments from registered C-arm images for long bone fracture reduction. In: IEEE-EMBS 2005: 27th annual international conference of the Engineering in Medicine and Biology Society, Engineering in Medicine and Biology Society, 1–4 Sept 2005, Shanghai, China. IEEE, Piscataway, NJ, pp 4361–4364 4. Hofstetter R, Slomczykowski M, Krettek C et al (2000) Computerassisted fluoroscopy-based reduction of femoral fractures and antetorsion correction. Comput Aided Surg 5(5):311–325. doi:10.3109/ 10929080009149849 5. Hüfner T, Pohlemann T, Tarte S et al (2001) Computer-assisted fracture reduction: novel method for analysis of accuracy. Comput Aided Surg 6(3):153–159. doi:10.1002/igs.1018 6. Joskowicz L, Milgrom C, Simkin A et al (1998) Fracas: a system for computer-aided image-guided long bone fracture surgery. Comput. Aided Surg. 3(6):271–288. doi:10.3109/10929089809148148 7. Nakajima Y, Tashiro T, Sugano N et al (2007) Fluoroscopic bone fragment tracking for surgical navigation in femur fracture reduction by incorporating optical tracking of hip joint rotation center. IEEE Trans Biomed Eng 54(9):1703–1706. doi:10.1109/TBME. 2007.900822 8. Oszwald M, Westphal R, Calafi A et al (2010) A standardized fracture reduction model for long bones-implication and evaluation in the femur. Technol Health Care 18(6):387–391 9. Ron O, Joskowicz L, Milgrom C et al (2002) Computer-based periaxial rotation measurement for aligning fractured femur fragments from ct: a feasibility study. Comput Aided Surg 7(6):332–341. doi:10.1002/igs.10056 10. Schmucki D, Gebhard F, Grützner PA et al (2004) Computer aided reduction and imaging. Injury 35(1):96–104 11. Tockus L, Joskowicz L, Simkin A et al (1998) Computeraided image-guided bone fracture surgery: modeling, visualization, and preoperative planning. In: Medical image computing and computer-assisted interventation, MICCAI’98, vol 1496, pp 29– 38. doi:10.1007/BFb0056185 12. Gösling T, Westphal R, Faülstich J et al (2006) Forces and torques during fracture reduction: intraoperative measurements in the femur. J Orthop Res 24(3):333–338. doi:10.1002/jor.20045 13. Graham AE, Xie SQ, Aw KC et al (2008) Bone–muscle interaction of the fractured femur. J Orthop Res 26(8):1159–1165. doi:10. 1002/jor.20611 14. Hung S, Lee M (2010) Functional assessment of a surgical robot for reduction of lower limb fractures. Int J Med Robot Comput Assist Surg 6(4):413–421 15. Koo TKK, Mak AFT (2007) A bone reposition device for execution of ct-based diaphyseal fracture reduction. J Biomech 40:S283 16. Matthews F, Neuhaus V, Schmucki D et al (2005) Passive pneumatic stabilization device for assisting in reduction of femoral shaft fractures. Eur J Trauma 31(6):568–574. doi:10.1007/ s00068-005-2047-3 17. Westphal R, Gösling T, Oszwald M et al (2008) Robot assisted fracture reduction. Experimental robotics, Springer tracts in advanced robotics, vol 39, pp 153–163. doi:10.1007/978-3-540-77457-0_15 18. Westphal R, Winkelbach S, Gösling T et al (2008) Telemanipulated long bone fracture reduction. In: Bozovic V (ed) Medical robotics. I-Tech Education and Publishing, Vienna, pp 507– 526 19. Westphal R, Winkelbach S, Wahl F et al (2009) Robot-assisted long bone fracture reduction. Int J Robot Res 28(10):1259–1278 20. Westphal R (2007) Sensor based surgical robotics. Contributions to robot assisted fracture reduction. Shaker, Aachen 21. Graham AE, Xie SQ, Aw KC et al (2006) Design of a parallel long bone fracture reduction robot with planning treatment tool. In: Proceedings of the 2006 IEEE/RSJ; international conference on intelligent robots and systems, pp 1255–1260

123

22. Graham AE, Xie SQ, Aw KC et al (2007) Robotic long bone fracture reduction. In: Bozovic V (ed) Medical robotics. I-Tech Education and Publishing, Vienna, pp 85–102 23. Joung S, Kamon H, Liao H et al (2008) A robot assisted hip fracture reduction with a navigation system. In: Medical image computing and computer-assisted intervention, MICCAI 2008. Springer, pp 501–508 24. Ye R, Chen Y (2009) Development of a six degree of freedom (DOF) hybrid robot for femur shaft fracture reduction. In: Proceedings of the 2008 IEEE, international conference on robotics and biomimetics, pp 306–311 25. Warisawa S, Ishizuka T, Mitsuishi M et al (2004) Development of a femur fracture reduction robot. In: IEEE international conference on robotics and automation, vol 4. IEEE, Piscataway, NJ, pp 3999– 4004 26. Füchtmeier B, Egersdoerfer S, Mai R et al (2004) Reduction of femoral shaft fractures in vitro by a new developed reduction robot system ‘RepoRobo’. Injury 35(1):113–119. doi:10.1016/j.injury. 2004.05.019 27. Mukherjee S, Rendsburg M, Xu WL (2005) Surgeon-instructed, image-guided and robot-assisted long bone fractures reduction. In: 1st international conference on sensing technology, pp 78–84 28. Seide K, Faschingbauer M, Wenzl ME et al (2004) A hexapod robot external fixator for computer assisted fracture reduction and deformity correction. IJMRCAS 01(01):64–69. doi:10.1581/ mrcas.2004.010101 29. Majidi Fakhr K, Kazemirad S, Farahmand F (2009) Robotic assisted reduction of femoral shaft fractures using stewart platform. Stud Health Technol Inform 142:177–179 30. Joung S, Shikh SS, Kobayashi E et al (2011) Musculoskeletal model of hip fracture for safety assurance of reduction path in robotassisted fracture reduction. In: Abu-Osman NA, Ting H (eds) 5th Kuala Lumpur international conference on biomedical engineering, 20–23 June 2011, Kuala Lumpur, Malaysia, vol 35. Springer, Berlin, pp 116–120 31. Kristen A, Culemann U, Fremd R et al (2008) Visualisierung von repositionspfaden. Unfallchirurg 111(6):395–402. doi:10.1007/ s00113-008-1429-5 32. Lorensen WE, Cline HE (1987) Marching cubes: a high resolution 3d surface construction algorithm. SIGGRAPH Comput. Graph. 21(4):163–169. doi:10.1145/37402.37422 33. Lohmann G (1998) Volumetric image analysis. Wiley, Teubner 34. Ohtake Y, Belyaev A, Seidel H (2004) Ridge-valley lines on meshes via implicit surface fitting. ACM Trans. Graph. 23(3):609. doi:10. 1145/1015706.1015768 35. Kühnel W (2008) Differentialgeometrie. Kurven-FlächenMannigfaltigkeiten, 4th edn. Vieweg, Wiesbaden 36. Carmo MPd (1976) Differential geometry of curves and surfaces. Prentice-Hall, Upper Saddle River 37. Monga O, Benayoun S (1995) Using partial derivatives of 3D images to extract typical surface features. Comput vis image underst 61(2):171–189. doi:10.1006/cviu.1995.1014 38. Therion J (1993) The extremal mesh and the understanding of 3D surfaces. Rapport de recherche, Institut National de Recherche en Informatique et en Automatique, vol 2149. INRIA, Le Chesnay 39. Lengagne R, Tarel J, Monga O (1996) From 2D images to 3D face geometry. In: Proceedings of the second international conference on automatic face and gesture recognition, Oct 14–16 1996, Killington. Vermont. IEEE Computer Society Press, Los Alamitos, pp 301–306 40. Kim S (2012) Extraction of ridge and valley lines from unorganized points. Multimed Tools Appl 63(1):265–279. doi:10.1007/ s11042-012-0999-y 41. Besl PJ, McKay ND, Schenker PS (1992) A method for registration 3-d shapes. IEEE Trans Pattern Anal Mach Intell 14(2):586–606. doi:10.1117/12.57955

Int J CARS 42. Guéziec A, Pennec X, Ayache N (1997) Medical image registration using geometric hashing. IEEE Comput Sci Eng 4(4):29–41. doi:10.1109/99.641607 43. Guéziec A, Ayache N (1991) Smoothing and matching of 3-D space curves. Rapports de recherche. Institut National de Recherche en Informatique et en Automatique, vol 1544. Institut National de Recherche en Informatique et en Automatique, Le Chesnay 44. Pajdla T, van Gool L (1995) Matching of 3-D curves using semidifferential invariants. In: 5th international conference on computer vision, IEEE Computer Society Press, pp 390–395 45. Pottmann H, Wallner J, Huang Q et al (2009) Integral invariants for robust geometry processing. Comput Aided Geom Des 26(1):37– 60. doi:10.1016/j.cagd.2008.01.002 46. Siciliano B, Khatib O (2008) Springer handbook of robotics. Springer, Berlin 47. Forsyth DA, Ponce J (2012) Computer vision. A modern approach, 2nd edn. Pearson, Boston 48. Winkelbach S (2006) Das 3d-puzzle-problem. Effiziente Methoden zum paarweisen Zusammensetzen von dreidimensionalen Fragmenten. Fortschritte in der Robotik, vol 10. Shaker, Aachen 49. Winkelbach S, Westphal R, Goesling T (2003) Pose estimation of cylindrical fragments for semi-automatic bone fracture reduction. In: Krell G, Michaelis B (eds) Pattern recognition, proceedings. Magdeburg, 10–12 Sept 2003, vol 2781. Springer, Berlin, pp 566– 573

50. Winkelbach S, Rilk M, Schönfelder C et al (2004) Fast random sample matching of 3D fragments. In: Rasmussen CE (ed) Pattern Recognition. 26th DAGM Symposium, Tübingen, Germany, August 30-September 1, 2004: proceedings, vol 3175. Springer, Berlin, New York, pp 129–136 51. Joskowicz L, Kronman A (2013) Automatic bone fracture reduction by fracture contact surface identification and registration. In: IEEE 10th international symposium on biomedical imaging, pp 246–249 52. Moghari MH, Abolmaesumi P (2008) Global registration of multiple bone fragments using statistical atlas models: feasibility experiments. In: Dumont G, Galiana H, Proncipe J et al (eds) Proceedings of the 30th annual international conference of the IEEE Engineering in Medicine and Biology Society. “Personalized Healthcare through Technology”; 20–24 Aug 2008, Vancouver Convention & Exhibition Centre, Vancouver, BC. IEEE Operations Center, Piscataway, NJ, pp 5374–5377 53. Albrecht T, Vetter T (2012) Automatic fracture reduction. In: Mesh processing in medical image analysis, vol 7599. Springer, pp 22–29

123