paulano-godino2016

Computer Methods and Programs in Biomedicine 141 (2017) 93–104

Contents lists available at ScienceDirect

Computer Methods and Programs in Biomedicine

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

Identification of fracture zones and its application in automatic bone fracture reduction

Félix PaulanoGodino ∗, Juan J. JiménezDelgado 1

Graphics and Geomatics Group, University of Jaén, Jaén, Spain

a r t i c l e i n f o

Article history: Received 29 March 2016 Revised 10 November 2016 Accepted 22 December 2016

Keywords: Alignment Ankle fractures Automatic bone fracture reduction Fracture zone Matching Preoperative planning

a b s t r a c t

Background and objective: The preoperative planning of bone fractures using information from CT scans increases the probability of obtaining satisfactory results, since specialists are provided with additional information before surgery. The reduction of complex bone fractures requires solving a 3D puzzle in order to place each fragment into its correct position. Computerassisted solutions may aid in this process by identifying the number of fragments and their location, by calculating the fracture zones or even by computing the correct position of each fragment. The main goal of this paper is the development of an automatic method to calculate contact zones between fragments and thus to ease the computation of bone fracture reduction.

Methods: In this paper, an automatic method to calculate the contact zone between two bone fragments is presented. In a previous step, bone fragments are segmented and labelled from CT images and a point cloud is generated for each bone fragment. The calculated contact zones enable the automatic reduction of complex fractures. To that end, an automatic method to match bone fragments in complex fractures is also presented.

Results: The proposed method has been successfully applied in the calculation of the contact zone of 4 different bones from the ankle area. The calculated fracture zones enabled the reduction of all the tested cases using the presented matching algorithm. The performed tests show that the reduction of these fractures using the proposed methods leaded to a small overlapping between fragments.

Conclusions: The presented method makes the application of puzzlesolving strategies easier, since it does not obtain the entire fracture zone but the contact area between each pair of fragments. Therefore, it is not necessary to find correspondences between fracture zones and fragments may be aligned two by two. The developed algorithms have been successfully applied in different fracture cases in the ankle area. The small overlapping error obtained in the performed tests demonstrates the absence of visual overlapping in the figures.

© 2016 Elsevier Ireland Ltd. All rights reserved.

1. Introduction

The development of computerassisted techniques increases the information available to specialists before surgery. This informa tion allows reducing surgery time and potential misinterpretations, with the consequent benefits to treatment and recovery time [1] . The reduction of a bone fracture requires properly identifying all the bone fragments, placing them to their original position and stabilizing the fracture. Computerassisted techniques facilitate the

∗Corresponding author. Present address: University of Jaén, Ed A3103, Jaén, Spain. Email addresses: fpaulano@ujaen.es (F. PaulanoGodino), juanjo@ujaen.es (J.J. JiménezDelgado).

1 Present address: University of Jaén, Ed A3142, Jaén, Spain.

process by enabling interaction with virtual models of fragments, by assisting in the planning of the surgery, by detecting the lack of tissue or by analyzing different configuration of fixation devices. The computerassisted preoperative planning of a fracture reduc tion may be divided into three steps: obtaining 3D models repre senting bone fragments, calculating the reduction of the fracture, and analyzing the obtained results [2] . The reduction of complex bone fractures requires solving a 3D puzzle in order to place each fragment into its proper position. In comminuted fractures, bone fragments usually have a manyto many relationship, that is, a fragment may share its fracture zone with more than one fragment. For this reason, additional process ing is necessary to determine which part of the fracture zone is shared with each fragment. Depending on the complexity of the

http://dx.doi.org/10.1016/j.cmpb.2016.12.014 01692607/© 2016 Elsevier Ireland Ltd. All rights reserved.

94 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

fracture, the surgery could last several hours. The previous iden tification of the bone fragments, their matching, and their correct positioning could help specialists to plan the surgery and thus to reduce the intervention time [1] . The main contribution of this paper is an automatic algorithm to calculate the contact zone between two fragments. The con tact zone may be defined as the area of the fracture zone that is shared between two fragments. Contrary to what happens with complete fracture zones, the calculated contact zones ease the use of puzzlesolving methods to compute the fracture reduction, be cause the correspondence between fracture zones is also implic itly calculated. The calculation of the contact zone between every pair of fragments is useful for the preoperative planning of a bone fracture reduction. On the one hand, the fracture reduction may be computed by matching and registering contact zones. On the other hand, contact zones provide additional information to better understand the fracture. The process of calculating contact zones can be summarized as follows. Firstly, bone tissue is segmented from a CT stack and each bone fragment is identified. In order to segment and label bone fragments from CT images, the approach proposed in [3] is used. The results of the segmentation are pro cessed to extract a point cloud for each fragment. Since the de veloped algorithm only requires the topological information of the contours, it is not necessary to generate 3D meshes representing bone fragments. Then the contact zone between fragments is cal culated using the method presented in this paper. The calculated contact zones can be matched and registered in order to reduce the fracture. For that purpose, a method to calculate the reduc tion of complex bone fractures using the obtained contact zones is also presented. The proposed procedure obtains good results in the case of moderate fracture displacement. If the fragments are too displaced from their original position, a previous coarse alignment of the fragments may be required. With that goal, the contralateral bone can be used as described in the literature [4,5] . Otherwise, in teractive methods can be utilized especially when the contralateral bone is not available [6,7] . The proposed method has been success fully applied in the reduction of different bone fractures. This paper is structured as follows. Section 2 reviews the cur rently proposed approaches to calculate the fracture zone and to reduce bone fractures. Section 3 describes the proposed methods for identifying contact zones and matching and aligning multiple bone fragments. After that, Section 4 tests the suitability of the proposed algorithms in the reduction of different bone fractures. The obtained results are discussed in Section 5 , extracting the main strengths and weaknesses of our approach. Finally, the paper is concluded in Section 6 .

2. Background

In the literature, different approaches have been proposed to compute the reduction of bone fractures. Some works propose to reduce the fracture by calculating, matching and registering the fracture zones. Alternatively, other authors propose interactive tools, or intend to reduce the fracture by registering the bone frag ments to a healthy template. In this section, the most relevant works are analyzed and classified. This classification helps to ex tract the main drawbacks of currently proposed methods.

2.1. Fracture zone calculation methods applied to virtual fracture reduction

Puzzlesolving methods require the calculation of the fracture zone of all the fragments involved in the fracture. Some differ ent methods have been proposed to calculate the fracture zone af ter the segmentation of bone fragments. Then the fracture can be virtually reduced by matching and aligning the calculated fracture

zones. Winkelbach et al. [8] take advantage of the specific shape of cylindrical bones. In order to identify vertices of the fractured area, they check the normal orientation of each vertex and compare it with the bone axis. This method does not work when fracture lines are almost parallel to the bone axis. Bone fragments are matched and registered by applying a modified Iterative Closest Point (ICP) algorithm to the fracture areas. Statisticalbased approaches have been proposed to identify fractured zones. Willis et al. [9] use a mixture model consisting of two Gaussian probability distributions to perform a binary classi fication that allows separating intact and fractured zones of each bone fragment. After classifying all points, the fractured surface is the largest continuous region of fractured surface points. Then the user has to interactively select fracture surface patches in pairs that coarsely correspond. Finally, bone fragments are automatically aligned using a modified ICP algorithm. Zhou et al. [10] present an extension of the previous method that improves fragment align ment in highly fragmented bone fractures. In order to separate fractured and intact surfaces, they use a twoclass Bayesian clas sifier based on the intensity values previously mapped on the sur face vertices. They upgrade the interactive method introduced in [9] by providing a userdirected search in order to match the frag ments. To perform the final alignment of the fracture surfaces, au thors also describe a new algorithm to improve the registration proposed by Willis et al. [9] . Other authors have proposed interactive methods to identify fracture surfaces in craniofacial fractures [11,12] . In these works, fracture contours are extracted interactively from segmented bone fragments. With that aim, specialists have to select points be longing to the fractured area and then a contour tracing algo rithm generates the rest of the points. Once the fracture contours are calculated, the 3D surface is generated by collating the con tours extracted from each slice. Chowdhury et al. [11] propose to match fracture surfaces using a Maximum Weight Graph Matching (MWGM) algorithm. To align the fracture surfaces, authors propose a modification of the ICP algorithm. A variation of this approach is presented in [12] . In order to match the fractured areas, they formulate a matrix score based on the appearance of mandibular fragments in the CT image sequence. Then they propose the Max imum Cardinality Minimum Weight (MCMW) algorithm for estab lishing the correspondence between each pair of fracture surfaces while registering them using the ICP algorithm. Curvature analysis has also been used to identify fractured sur faces. Okada et al. [4] present a curvaturebased procedure to ob tain fracture lines in each slice. For that purpose, they generate a 3D curvature image from the CT stack. Then interactive line tracking software allows extracting the fracture lines from the gen erated 3D curvature image. In order to reduce proximal femur frac tures, authors propose to approximate bone fragments to their cor rect position using the contralateral bone as template and then finetune the result by registering the calculated fracture lines. That registration is performed using the ICP algorithm. Fürnstahl et al. [5] use a normalbased filter to identify points candidate to belong to the fracture surface. Then connected component anal ysis is applied in order to remove outliers. User interaction may be required to remove incorrect parts not belonging to the frac ture surface. Then authors use the contralateral bone to perform a coarse initial alignment of the bone fragments. Finally, frac ture lines are registered using both a pairwise registration and a multipiece alignment. Kronman and Joskowicz [13] identify frac ture surfaces using intensity and curvature filters. Then outliers are removed from the estimated contact surface with a connectivity test. The calculated fracture surfaces enable the reduction of sim ple fractures by applying ICP registration to the fracture surfaces. Buschbaum et al. [14] identify fracture lines using surface curva tures. Then the target position is calculated by matching fracture

F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104 95

lines. Finally, authors propose to use referencecoordinate system in order to compute reduction parameters. Currently proposed puzzlesolving methods calculate the com plete fracture zone of each bone fragment. Moreover, some of these methods require user interaction. The method presented in this pa per automatically calculates the contact zone between each pair of fragments. Therefore, no additional calculations are needed to ob tain correspondences between fragments and thus matching meth ods can be easily applied in order to reduce the fracture.

2.2. Other approaches to virtually reduce bone fractures

Not all the proposed computerassisted methods proposed to reduce bone fractures are based on the calculation of the frac ture zone. Once virtual bone fragment models have been gener ated, some authors proposed interactive tools to place them into their correct position. Cimerman et al. [15] present an application to interactively reduce pelvic and acetabular fractures. Acetabular fractures can also be reduced by using the software presented in [7] . Other authors have developed applications to manipulate bone fragments reconstructed from CT scans [16,17] . Haptic devices have been used for the interactive reduction of acetabular fractures [6] . The virtual reduction of a bone fracture is a complex procedure and thus interactive methods are time consuming. Furthermore, it is difficult to obtain accurate results if the final alignment is man ually performed. In order to make the reduction process as automatic as possi ble, some authors propose to perform a registration between the bone fragments and an intact template. These approaches do not require a matching procedure since the templates determine the position and orientation of each bone fragment. With the aim of defining the intact template, Maubleu et al. [18] propose to mir ror the healthy side of the bone. By using this procedure, authors simulate the reduction of a fractured zygomatic bone. Other au thors use a statistical anatomical atlas as a template [19] . The atlas is generated using principal component analysis (PCA) in a bone population generated from CT scans. Gong et al. [20] also utilize an anatomical atlas to register fragments of the distal radius. In the first place, each segmented volume is properly positioned on the atlas. After that, the atlas is deformed to maximize the over lap area between its 2D projection and the individual fragments in the segmented regions. Albrecht and Vetter [21] use the ICP algorithm to make a rigid alignment of the fragments and then they adapt the statistical shape model to the bone fragments. An other possibility is to use artificially generated bones as templates [22] . In this work, authors generated identical tibias from blocks of highdensity polyurethane foam. These tibias were scanned before and after being fractured in order to use them as a template. Kato [23] presents a framework for correspondenceless registration of 3D shapes. Using that framework and a template, they place each fragment into its correct position by formulating the problem as an affine puzzle. Templatebased works do not require user inter action but they are limited to the fracture cases defined by the templates; hence it is laborious to extend them to deal with dif ferent bones and fracture types. Moreover, the obtained results de pend on the quality of the template and not always it is possible to find a template that enables the registration. Recent studies have shown that the contralateral bone should be used with caution as a template, since differences in shape may exist between symmet rical bones [24] .

3. Methods

The presented methods can be divided into three main stages. Firstly, bone fragments are segmented from CT scans. The second stage consist of the calculation of contact zones between bone

fragments. For this purpose, a point cloud representing each bone fragment is extracted from the segmentation results. After that, a discretized sweep is carried out in order to obtain a set of points candidate to belong to the contact area. Since not all the candidate points belong to the fracture zone, the set of points is then filtered to get the actual contact zone. In the last stage, the calculated con tact zones are matched and aligned in order to reduce the fracture. During this stage, contact zones are recalculated after aligning each pair of bone fragments.

3.1. Bone fragment identification

The identification of bone tissue from CT images is a complex task; hence it is difficult to find a solution that works in all cases. In a bone, two very distinct zones can be distinguished: corti cal and trabecular tissue. Cortical tissue is very dense and can be found in the outer part of the bone. Trabecular tissue is mainly in the inner part of the bone, is more heterogeneous and has less in tensity in a CT image. Nevertheless, the cortical zone is very thin near joints and even disappears in the area closest to the join. Thus the transition of the intensity values near the joints gener ally appears to be fuzzy and some areas within the bone may have similar intensity than the soft tissue surrounding the bone. Dur ing the segmentation, this may cause incomplete segmentation or overgrowing. Fractured bone is even more difficult to be identified [25] , because additional constraints need to be considered. Bone fragments may have arbitrary shape and can belong to any bone in a nearby area because of the fracture. Therefore, all the fragments need to be labelled during or after the segmentation process. In addition, bone fragments are not completely surrounded by corti cal tissue due to the fracture. This is even more problematic near joints because of the aforementioned special distribution of bone tissue in this zone. Finally, the resolution of the CT image and the proximity between fragments may cause that different fragments appears as only one in the image. This may require additional pro cessing to split erroneously joined bone fragments. In order to segment and label bone fragments from CT images, the method proposed in [3] is used ( Fig. 1 ). This method is based on 2D region growing and has been successfully tested in the iden tification of bone fragments in different types of fractures. The user only has to identify bone fragments in the first slice in which each one appears. To that end, a seed must be placed inside the region of each target bone fragment. These seeds enable the execution of the segmentation algorithm and the labeling of the bone frag ments. Then 2D region growing is executed for each seed and after that all the seeds are propagated through the image stack. Unlike other 3Dbased approaches, this method is able to detect and solve overgrowing cases in the first slice they appear, since 2D region growing is applied in each slice. In order to solve these cases, 2D isolated region growing is proposed. This method allows splitting erroneously joined bone fragments during the segmentation pro cess. The user only has to place an additional seed in the fracture area ( Fig. 2 ). Then the algorithm calculates an intensity thresh old value that enables the separation of the fragments. Sometimes, there is no threshold value that separates the wrongly joined frag ments. For example, the case of two fragments connected by their cortical areas. For these cases, the method uses the information from the previous slice to estimate the shape of the fragment. The segmentation algorithm is fully described in [3] .

3.2. Calculation of the contact zone between fragments

With the aim of easing the application of puzzlesolving meth ods to reduce the fracture, the presented method does not cal culate all the fracture zones of a fragment. Instead, it calculates

96 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

Fig. 1. Point clouds representing bone fragments identified by the method described in [3] .

Fig. 2. Left, two bone fragments are wrongly joined as a result of applying 2D re gion growing. Right, result of the application of the isolated region growing based method.

the contact zone between each pair of fragments; hence the corre spondences between fragments are implicitly calculated. For simplicity, let us consider the calculation of the contact zone between two fragments. Given two bone fragments A and B , the method presented in this section calculates the area of A that should be joined with the area of B to reduce the fracture and vice versa, namely the contact zones CZ A ^{and CZ} B ^{between the two frag} ments. Section 3.3 shows how to use the calculated contact zones in order to compute the reduction of complex fractures.

3.2.1. Generation of point clouds representing bone fragments The input of the contact zone calculation algorithm is a pair of point clouds representing each bone fragment ( Fig. 1 ). These point clouds are calculated from the CT image stack resulting of the identification of bone fragments explained in Section 3.1 . Each slice of the stack contains regions of pixels representing the seg mented bone tissue. For each pixel in each image, a point is gener ated if the intensity value associated with the pixel is greater than 0. Therefore, the generated point clouds represent both cortical and trabecular tissue. Fig. 1 shows a cross sectional view of a generated bone fragment. In addition to the coordinates of the position of the point in 3D space, the data structure of each point stores the in tensity value of the corresponding voxel in the CT image stack and its associated normal vector. In a first approach, we tried to estimate normals in 3D by using the algorithm implemented in the Point Cloud Library (PCL) [26] .

In order to estimate the normal vector of a point, this algorithm uses the points that are closer than a specified Euclidean distance to it. Then the algorithm estimates the normal of the plane tan gent to the surface defined by the nearest points. However, the ob tained results were not accurate, since the distance between slices in the original CT stack was at least twice the size of pixels. As an alternative, normals are estimated in 2D for each point using topological information from contours. Thus, the normal vectors are only estimated at the contour points of the fragment. In a first stage, external contours are extracted from the segmented regions by applying Marching Squares (MS) [27] . Due to the fact that seg mented regions have holes, MS also obtains internal contours. In order to remove those internal contours, the pointinpolygon al gorithm by Feito and Torres [28] is used. Then normal vectors of the contour points are estimated. For each three consecutive points P i −1 , P i ^{and P} i +1 ^{in a contour, the normal associated to P} i ^{is esti} mated as the sum of the normals of the two segments formed by P i −1 ^−P i ^{and P} i ^−P i +1 . The third component of each normal is de fined by the plane of the slice. Despite of this estimation, the per formed tests showed that normals are accurate enough to estimate the surface curvature. Fig. 3 shows the estimated normals of two different fragments.

3.2.2. Generation of the candidate points The proposed algorithm performs a sweep to calculate points candidate for belonging to the contact zone between the two frag ments. The set of candidate points contains most of the points be longing to the contact zone but not all the candidate points be long to the contact zone. Our method carries out a sweep with the purpose of calculating the candidate points of each fragment. Therefore, the first step of the algorithm consists of the identifica tion of the sweep direction. For that purpose, the centroid of each point cloud must be calculated. Both centroids will be named C A and C B ^{respectively. The segment that joins the two centroids de} fines the direction of the sweep ( Fig. 4 , left). Since a point cloud is used to represent each bone fragment, the space must be dis cretized before performing the sweep. Candidate points are those that are reached first by the sweep in each of the spatial divisions. The algorithm is explained in detail in the following paragraphs. In order to perform the sweep, the volume of each point cloud is discretized using an Oriented Bounding Box (OBB) [29] for each point cloud. In order to simplify, only the process to calculate the contact zone CZ A ^{is explained. Once CZ} A ^{is obtained, the process} needs to be entirely repeated changing the sweep direction in order to obtain CZ B . The only restriction of the OBBs is that one of their associated direction vectors is determined by the segment that joins the centroids of both point clouds. Thus the orientation

F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104 97

Fig. 3. Front and top view of the estimated normals of point clouds representing two different bone fragments.

of the OBB of fragment A will be defined by the basis ( ⃗ x , ⃗ y , ⃗ z ) , where ⃗ x is the normalized result of C A ^−C B ^{and ⃗} x , ⃗ y and ⃗ z form an orthonormal basis. ( Fig. 4 , left). After being calculated, each OBB is homogeneously subdivided to build a grid, so that each voxel of

the grid stores all the points that are located inside it. The size of the grid voxels is defined by the distance between slices in the original stack. This assures that a small number of points – no more than two or three – are classified in each voxel of the grid. The data structure of the grid only stores the size of its OBB and the number of subdivisions in each dimension. Moreover, the re lationship between each point and its corresponding position on the grid is also saved. Therefore, there is no memory problems, since the maximum number of classified points is defined by the number of segmented pixels in the CT stack. Moreover, the grid structure is temporal and it is deleted after the calculation of the candidate points; hence there are only two grids concurrently in memory. Afterwards, a discretized sweep is performed using the grid. The direction of the sweep is from C B ^{to C} A . For each fixed value of x and y, the content of the grid voxels V xyz k ^{of fragment A, where k} goes from 0 to l, is studied in each step of the sweep. The sweep ing begins with all the voxels whose last coordinate is z 0 . If the voxel V xyz k ^{does not contain any point, it is discarded and the voxel} V xyz k +1 ^{is studied. On the contrary, if the voxel V} xyz k ^{contains at} least a point, all the points classified in that voxel are considered as contact zone candidates and the remaining voxels with these x and y coordinates are not checked. This process is illustrated in Fig. 4 , right. After finishing the sweep, a set of candidate points is obtained. Due to the size of subdivision chosen for the grid, the contact zone obtained practically represents a surface. This will allow the application of a registration algorithm to join the frag ments.

3.2.3. Filtering of the candidate points The set of candidate points needs to be filtered in order to ob tain the actual contact zone of each fragment ( Fig. 5 ). The revi sion performed in Section 2 showed that curvaturebased methods

Fig. 4. Left, OBBs and centroids associated with two fragments. Sweep direction used to calculate the fracture zone of fragment A. Right, 3D representation of the sweeping procedure. Checked voxels that contain no points are displayed in grey. Checked voxels that contain points candidate to belong to the contact area are displayed in orange. White cells are not checked.

98 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

Fig. 5. Left, set of candidate points of fragments A and B before applying any filter. Centre, fracture zone after the application of the distancebased filter. Right, frac ture zones after performing the registration with the aim of joining the two bone fragments.

are wellknown and frequently used to identify the fracture zone of bone fragments. The small number of clinical cases of fractures at our disposal and their variety make impossible the utilization of statisticsbased methods. Consequently, two different parame ters are considered for filtering contact zones: the distance from each point to the opposite fragment and the estimated curvature at each point. The first parameter is based on the minimum distance from each point to the other fragment. The complexity of computing the minimum distance would be O ( n 2 ), because each point should be tested with all the points of the opposite fragment. However, this value is calculated as the minimal distance between the point

and the candidate points of the opposite fragment, reducing com plexity. Since the contact zone must be close to the fragment to be joined, points whose minimum distance to the other fragment is higher than a predefined threshold must be discarded. This threshold depends on each individual case and it is established by experience. The second parameter is the estimated curvature at each point ( Fig. 6 ). Unlike the distancebased parameter, the estimation of the curvature requires the normal vector at each candidate point. Therefore, this filter only considers the contour points of each frag ment. The implementation available in PCL is used [26] in order to estimate the curvature at each contour point. This implementation projects the normals of the closest points onto the tangent plane in each point and then uses Principal Component Analysis (PCA) in order to estimate the principal curvature. It is expected that con tact zone points have a high curvature because the area is very irregular due to the fracture and it is mostly composed by tr abec ular tissue. The application of these filters allows discarding points result ing from the sweep algorithm that do not belong to the contact zone. The algorithm proposed for calculating the contact zone be tween two fragments is depicted in Fig. 7 . The utilization of one parameter or the other for filtering could depend on the bone and the fracture type. In Section 4 both filters will be evaluated in the reduction of different bone fractures.

3.3. Matching and alignment of multiple bone fragments

The reduction of a bone fracture requires placing all the bone fragments at their original positions. In simple fractures, only two fragments are generated; hence the problem is reduced to align their fracture zones. The contact zones calculated by the method proposed in this paper enable the application of the ICP algorithm

Fig. 6. Estimated maximum principal curvature of two different point clouds representing bone fragments. High curvature values are shown in blue and low curvature values are shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104 99

Fig. 7. Flowchart for the automatic calculation of the contact zone between two fragments.

[30] for registering the points belonging to the contact zones. The transformation matrix obtained by executing the ICP algorithm is applied to one of the fragments and as a result the fracture is re duced ( Fig. 5 , right). Fractures with more than two fragments need an additional matching procedure to establish the correspondence between frag ments. Since the method described in Section 3.2 calculates the contact zone between each pair of fragments, bone fragments can be incrementally matched in pairs. In order to establish the order in which each bone fragment pair is joined in the fracture reduc tion process, a matching score is calculated for each pair of frag ments. Thereby, the pair of fragments with the highest score has priority. As the goal is to ensure that the best joining is carried out in each cycle of the loop, the fitness measure obtained by the ICP algorithm is used as score. This fitness measure is calculated as the average of the squared minimum distances from the points of one contact surface to the other. Therefore, the first step consists of the identification of the con tact area and the calculation of the fitness measure between each pair of fragments. Then the two fragments that obtain the best score are joined. To that end, the transformation matrix obtained by the ICP algorithm is used as in the case of simple fractures. After that, the contact zones between the resulting fragment and

Fig. 8. Flowchart for reducing complex bone fractures using the proposed algorithm to calculate contact zones.

each one of the fragments that remain to be repositioned are cal culated. The entire process is then repeated until all the fragments have been matched and joined and thus the bone fracture has been reduced. This algorithm is schematized in Fig. 8 .

4. Results

The algorithms described in Section 3 have been successfully applied in the reduction of different bone fractures in the an kle area. Firstly, distancebased and curvaturebased filters have been tested with the aim of determining which filter produces bet ter results in the calculation of the contact zone for the available datasets. The threshold values for both filters are variable and they have been manually selected by experience in each case to avoid compromising the comparison. Subsequently, the performance of the proposed fracture reduction method has been measured. To that aim, the overlapping between fragments has been calculated and the runtime of the different parts of the method has been measured. The datasets contain three different fractures in the ankle area ( Fig. 9 ). The method proposed in [3] was able to segment and label bone fragments in all tested cases. The first CT stack con tains a complex distal fracture of tibia and fibula ( Fig. 9 , A), and it is formed by 173 slices with dimensions 159 × 159 mm and with 0.625 mm spacing. The fracture generated five bone frag ments. The second case is an irregular extraarticular fracture of talus that generated two fragments ( Fig. 9 , B). The CT stack con tains 156 slices. The dimensions of each slice are 150 × 150mm and the separation between slices is 0.625 mm. The third CT stack represents an extraarticular fracture of calcaneus ( Fig. 9 , C) and contains 205 slices. The fracture produced two fragments. The sep aration between slices is 0.625 mm and the dimensions of each

100 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

A C B

Fig. 9. Front and top view of the different fracture cases used for experiments. (A) Complex fracture of fibula and tibia, (B) simple fracture of talus and (C) calcaneus. All the bone fragments were identified using the method described in Section 3.1 .

Table 1 Parameters of the CT stacks used as input for experiments.

# slices Size of Resolution Spacing between each slice (mm) of slices slices (mm)

Fibula and tibia 173 159 ×159 512x512 0.625 Talus 156 150 ×150 512x512 0.625 Calcaneus 205 173 ×173 512x512 0.625

slice are 173 × 173 mm. The resolution of the images of the three stacks is 512 × 512. Table 1 summarizes the parameters of all the CT stacks used as input for experiments.

4.1. Testing different parameter configurations

The two proposed filters were tested for obtaining the con tact zone between different bone fragments. The first test uses the curvature filter and the second test uses the distancebased filter. From this point, they will also be called curvaturebased and distancebased method respectively. Curvature and distance thresholds were manually established for each particular case in order to obtain the best result and thus to avoid compromising the comparison between both methods. In both cases, an additional ra dius outlier removal filter was applied in order to remove isolated points. For the first fracture, the curvaturebased method did not per form well because some points in the joint area have similar curva ture values than the points in the fracture zone ( Fig. 10 , A, bottom). By using the distancebased method, the contact zones were cor rectly calculated ( Fig. 10 , A, top). Fig. 11 shows that the calculated contact zones enabled a proper reduction of the fracture of fibula using the distancebased method, but not using the curvature based method. The fracture of tibia required applying the match ing algorithm in order to establish the two bone fragments to be joined in the first place. After joining the first pair of fragments, the contact zone of the resulting fragment was calculated with re

spect to the remaining fragments. Finally, the fracture was com pletely reduced by registering these two new contact zones. Due to the irregular shape of the bone in the talus fracture, the curvature filter performed even worse than in the previous case. Many points in the outer part of the bone have an estimated curvature similar to the points in the fracture zone; hence the curvaturebased filter obtained a scattered point cloud ( Fig. 10 , B, bottom). Nevertheless, the distancebased filter was able to auto matically compute contact zones ( Fig. 10 , B, top). Then the fracture was reduced by registering them. The irregular shape of the bone in the calcaneus fracture caused that the curvature method also obtained a scattered point cloud ( Fig. 10 , C, bottom). However, the distancebased filter was able to deal with the calculation of the contact zone automatically ( Fig. 10 , C, top). Then the registration algorithm completed the fracture re duction.

4.2. Performance of the proposed method

The results of the proposed method using the distancebased filter have been quantitatively measured. The quality of curvature based filter has not been measured, since it obtained bad re sults in the qualitative tests. The ideal would be to compare the obtained results to the ground truth, but this is not available for our clinical cases. Alternatively, the overlapping between frag ments has been considered, since it is a test parameter commonly used in the literature [4,5] . After performing the fracture reduction

F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104 101

A C B

A C B

Fig. 10. Fracture zones obtained by applying the distancebased filter (top) and the curvaturebased filter (bottom). (A) Complex fracture of fibula and tibia, (B) simple fracture of talus and (C) calcaneus.

Fig. 11. Results obtained for the fracture of fibula using the distancebased filter (left) and the curvaturebased filter (right), including the calculated fracture zones.

process, bone fragments should be in touch. Fig. 12 shows that there was no visible gap between fragments after the fracture re duction process. Nonetheless, tests have been performed to quan titatively measure that bone fragments were in touch but did not overlap. Since bone fragments were represented by point clouds, they were discretized before calculating the overlapping between fragments. For that purpose, first the AxisAligned min imum Bounding Box (AABB) [29] of the resulting point cloud was computed. Then the AABB was subdivided in voxels of the same size as the original voxels of the CT image stack. With the goal of enabling subsequent overlap calculations, each voxel referenced points that were located inside and each point stored the bone fragment to which it belonged.

In order to compute the overlapping, all the AABB voxels were processed. Touching voxels were defined as those containing points belonging to two different fragments. Then all the touching voxels were candidates to be overlapping voxels. The calculation of over lapping voxels between two fragments required studying the 26 neighbors of each touching voxel. If there is overlapping, touching voxels are surrounded by other touching voxels. For this calcula tion, empty voxels were not considered. Therefore, touching vox els were also overlapping voxels if the number of their touching neighbor voxels was greater than 2/3 of the number of their non empty neighbor voxels. In the case of the complex fracture of tibia, the process was repeated for each pair of fragments and the ob tained overlapping voxels were accumulated. The overlapping error

102 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

A C B

Fig. 12. Fractures reduced by applying the distancebased method. (A) Complex fracture of fibula and tibia, (B) simple fracture of talus and (C) calcaneus.

Table 2 Efficacy of the distancebased fracture reduction process. The overlapping error is defined as the ratio between overlapping and touching voxels.

#Touching #Overlapping Overlapping/Touching cells cells ratio (Overlapping error)

Fibula 218 8 0.0367 Tibia 1683 148 0.0879 Talus 210 12 0.0571 Calcaneus 105 4 0.0381

was defined as the ratio between overlapping voxels and touching voxels. Table 2 shows the results. The obtained overlapping error for the three tested cases was very small. Although the proposed method was able to compute the frac ture reduction without user interaction in the three tested cases, it may require a previous coarse alignment of the bone fragments in some cases. In order to carried out the previous alignment, our developed interactive application [17] can be used. Since it is dif ficult to get accuracy when manually aligning bone fragments, our automatic methods should be robust to small translations and ro tations. In order to assess the robustness of our method, we ran domly translated bone fragments up to 3 mm and rotated them up to 5 degrees in each dimension before applying it. We con ducted 15 experiments for each tested case, discarding the exper iments in which bone fragments were colliding in their initial po sitions. The computed final position of the previously translated and rotated bone fragment was compared to its computed final position without being previously modified. The translation error was calculated as the Euclidean distance between the centroid of the reduced bone fragment in its computed final position being and without being previously translated and rotated. In the same way, the rotational error was calculated as the average difference around the fragments three 2nd moment vectors. The obtained av erage translation was 1.0892 ± 0.546 mm and the rotation error was 3.1242 ± 1.3775 degrees, proving that the proposed method is quite robust to small variations in the initial position of the fragments.

Table 3 Translation and rotation errors obtained by the proposed method in comparison to virtual reductions carried out by a specialist.

Translation error (mm) Rotation error (degrees)

Fibula 1.52212 4.3658 Tibia 1.61343 2.4283 Talus 1.64876 4.192 Calcaneus 2.39704 2.0048

With the aim of complementing the assessment of the accuracy of the proposed method, the obtained reduced fracture could be compared to an intact template. However, the contralateral bone and the bone before being fractured are not available for the tested clinical cases. Alternatively, a specialist virtually reduced the 4 tested cases using our interactive application [17] in order to al low the assessing of our method by comparing it to the virtual reduction. This procedure provides a measure of how similar are the results obtained to the specialist criterion. As in the previous experiment, the translation error was defined as the Euclidean dis tance between the centroid of the fragment in the automatically computed and the virtually reduced position. Similarly, the rota tional error was calculated as the average difference around the fragments three 2nd moment vectors. Translation error was above 2.5 mm and rotation error was smaller than 5 degrees on average in the four tested cases ( Table 3 ). Additionally, the efficiency of each step of the algorithm was measured ( Table 4 ). The tests were run on a PC with one first gen eration Intel Core i7 2.8Gz and 4GB RAM. The contact zone was calculated using the distancebased method. Once contact zones were calculated, the time to register them was also measured. The table includes all the times each stage was run during the fracture reduction process. In the case of simple fractures, each stage needs to be executed once. For the tibia fracture, contact zones were cal culated 4 times and the registration process was executed twice. Table 4 shows that the fracture reduction process requires a few seconds for simple cases and over a minute for the more complex

F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104 103

Table 4 Efficiency of the distancebased fracture reduction process. The tests were run on a PC with one first generation Intel Core i7 2.8Gz and 4GB RAM. Runtime is given in seconds.

Calculate #Fracture Register #Fracture Fracture #points fracture zone fracture zone reduction zone (s) calculations zones (s) registrations (Total time) (s)

Fibula 3 .15 1 0.04 1 3 .21 125599 Tibia 18 .70 4 0.06 2 74 .87 401416 Talus 8 .37 1 0.09 1 8 .49 519812 Calcaneus 11 .39 1 0.08 1 11 .51 346087

case. These runtimes should not be an obstacle to its utilization in the preoperative planning of a bone fracture.

5. Discussion

This paper has presented a method to automatically calculate the contact zone between each two fragments. The calculated frac ture zones enable the reduction of complex bone fractures. The main advantage of our method is that it calculates the contact zone between each pair of fragments instead of the entire fracture zone, thus it properly works without modification if any fragment is missing as long as the fracture can be reduced. Unlike other meth ods proposed in the literature [4,5,14] , our method avoids the cal culation of correspondences between fragments, which is a com plex and timedemanding procedure. Moreover, the method pro posed in this paper makes use of both cortical and trabecular tis sue in order to join bone fragments. This causes that the contact areas between fragments are greater, allowing a better alignment of the bone fragments. Sometimes, trabecular tissue cannot be seg mented from CT images due to its low intensity or even by the parameters of the CT scanner. Other times, trabecular tissue can be deformed because of the fracture. In order to be able to prop erly compute the fracture reduction in these cases, the proposed method should be adapted to work with only cortical tissue. Regarding the final alignment, our approach and most of the currently proposed works [2] use an adaptation of the ICP algo rithm. Although it was not necessary in the conducted experi ments, the method proposed to calculate the contact zone in this paper may require a previous coarse alignment in the case of displaced fractures, especially if bone fragments are significantly moved or rotated with respect to their original position. This oc curs particularly in the case of crush injuries. However, the tests conducted demonstrated that our method is robust to small dis placements in the initial position of the bone fragments, thus the previous alignment does not need to be too accurate. In order to carry out a previous alignment of the fragments, a templatebased rough alignment is usually executed before applying the fracture reduction method [4,5] . The use of templates is limited, since it is not easy to get them in clinical cases. As an alternative, we devel oped an interactive application [17] to carry out a proper previous alignment of bone fragments when needed. Our clinical cases contain 3 different fracture cases in which 4 different bones are involved. This small number of cases is due to the fact that it is difficult to obtain CT scans of fracture cases that can be reduced without the use of osteoimplants. This is a recur rent problem in the literature. [4,5,10] . On the other hand, there is no standard criterion for evaluating the accuracy of the fracture reduction composition. In the literature, authors have proposed dif ferent methods to measure the quality of the obtained fracture re duction; hence the different proposed methods cannot be easily compared. In most cases authors used artificially generated bones to test their proposed techniques [10,20,22] . In these cases, the re sults obtained cannot be compared, since they used different bones to perform the tests: corpse bones and artificial bones of different

materials. The special features of long bones have also been used to check the accuracy of the reduction [31] . These features enable the comparison between different methods, but this is restricted to long bone fractures. Registrationbased methods can use the Mean Squared Error (MSE) to test the accuracy of the final alignment of the fracture surfaces [12,32] . This measurement can only eval uate this step of the process, it is only applicable in methods that use registration techniques, and it allows comparing the registra tion of the same surfaces. With the aim of making comparatives easier, the fracture reduction method presented in this paper has been evaluated using parameters commonly used in the literature to evaluate clinical cases [4,5] . Specifically, overlapping and con tact surfaces were calculated. In addition, the results obtained by our method have been compared to virtual interactive fracture re ductions carried out by a specialist, providing a measure of how similar are our results to the specialist criterion. The conducted tests show that contact zones calculated by the distancebased method enable the application of the proposed matching method to reduce different fractures in the ankle area. On the other hand, the curvature does not seem to be a good parameter to filter fracture points of bone fragments using our models. However, this parameter has been successfully applied in the literature to filter points belonging to the fracture zone in model representing only cortical tissue [4,5,14] . This may be be cause these methods only consider cortical tissue for the calcula tion of the fracture zone. In order to obtain acceptable results in our test cases, the curvaturebased method would need user inter action to discard points outside the fracture zone.

6. Conclusions

In this paper, an algorithm to automatically calculate the con tact zone between bone fragments is presented. The calculated contact zones enable the automatic reduction of complex fractures. To that end, an automatic method to match bone fragments in complex fractures has also been proposed. The approach presented in this paper has been successfully ap plied in different fracture cases in the ankle area. The obtained re sults show that fracture reductions have been successfully com puted for all the tested clinical cases, since there is no visible overlapping or gaps between fragments. The small overlapping er ror demonstrates the absence of visual overlapping in the figures. Moreover, conducted tests demonstrated that our method is ro bust to small displacements in the initial position of the bone frag ments. Additionally, the results obtained by the proposed method have been compared to virtual fracture reductions carried out by a specialist by using an interactive tool. These tests have provided a measure of how approximate are our results to the specialist crite rion. The calculation of the contact zone only takes a few seconds for the more complex cases. Nevertheless, this calculation could be easily parallelized in multicore architectures due to the fact that the algorithm is based on the individual processing of each point of the bone fragments. In this way, overall time could be consid erably reduced. Despite these results, a clinical study should be

104 F. PaulanoGodino, J.J. JiménezDelgado / Computer Methods and Programs in Biomedicine 141 (2017) 93–104

carried out to determine if the developed method would lead to an improvement in the preoperative planning of bone fracture re duction. In the future, new methods will be developed with the aim au tomatically select proper threshold values for the distance and cur vature filters. These methods would allow a full automation of the methods presented in this paper. Although we developed interac tive tools to carry out a previous rough alignment of bone frag ments, the ideal would be to develop new templatefree methods with the goal of performing an initial alignment that enables the use of the proposed method for reducing fractures containing very displaced bone fragments. Since the proposed method is not based on intensity, it could be adapted to only work with cortical tissue in the future. This modification would allow the method dealing with missing and deformed trabecular tissue. One of the main dif ficulties of current developments in the literature is the obtaining of an extensive and diverse set of clinical cases with which to eval uate the proposed methods. A possible solution would be to de velop fracture simulation methods that could provide on demand diverse realistic fracture cases of different bones from 3D models of healthy bones.

Acknowledgments

This work has been partially supported by the Ministerio de Economía y Competitividad and the European Union (via ERDF funds) through the research project DPI201565123R.

References

[1] N. Sugano, Computerassisted orthopedic surgery, J. Orthop. Sci. 8 (3) (2003) 4 42–4 48, doi: 10.1007/s10776 002 0623 6 .

[2] J.J. Jiménez, F. Paulano, R. Pulido, J. Jiménez, Computer assisted preoperative planning of bone fracture reduction: simulation techniques and new trends, Med. Image Anal. 30 (2016) 30–45, doi: 10.1016/j.media.2015.12.005 .

[3] F. Paulano, J.J. Jiménez, R. Pulido, 3D segmentation and labeling of fractured bone from CT images, Vis. Comput. 30 (68) (2014) 939–948, doi: 10.1007/ s00371 014 0963 0 .

[4] T. Okada, Y. Iwasaki, T. Koyama, N. Sugano, Y.W. Chen, K. Yonenobu, Y. Sato, Computerassisted preoperative planning for reduction of proximal femoral fracture using 3DCT data., IEEE Trans. Biomed. Eng. 56 (3) (2009) 749–759, doi: 10.1109/TBME.20 08.20 05970 .

[5] P. Fürnstahl, G. Székely, C. Gerber, J. Hodler, J.G. Snedeker, M. Harders, Com puter assisted reconstruction of complex proximal humerus fractures for pre operative planning., Med. Image Anal. 16 (3) (2012) 704–720, doi: 10.1016/j. media.2010.07.012 .

[6] J. Fornaro, M. Keel, M. Harders, B. Marincek, G. Székely, T. Frauenfelder, An interactive surgical planning tool for acetabular fractures: initial results., J. Or thop. Surg. Res. 5 (50) (2010) 1–8, doi: 10.1186/1749 799X 5 50 .

[7] Y. Hu, H. Li, G. Qiao, H. Liu, A. Ji, F. Ye, Computerassisted virtual surgical pro cedure for acetabular fractures based on real CT data, Injury 42 (10) (2011) 1121–1124, doi: 10.1016/j.injury.2011.01.014 .

[8] S. Winkelbach, R. Westphal, T. Goesling, Pose estimation of cylindrical frag ments for semiautomatic bone fracture reduction, in: B. Michaelis, G. Krell (Eds.), Pattern Recognition, LNCS 2781, Lecture Notes in Computer Science, vol. 2781, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 566–573, doi: 10.1007/978 3 540 45243 0 _ 72 .

[9] A. Willis, D. Anderson, T. Thomas, T. Brown, J.L. Marsh, 3D reconstruction of highly fragmented bone fractures, in: Proceedings of the SPIE. Medical Imaging 2007: Image Processing, 2007, p. 65121, doi: 10.1117/12.708683 .

[10] B. Zhou, A. Willis, Y. Sui, D. Anderson, T.P. Thomas, T. Brown, Improving inter fragmentary alignment for virtual 3d reconstruction of highly fragmented bone fractures, in: SPIE Medical Imaging, 2009, pp. 7259341–7259349, doi: 10.1117/ 12.810967 .

[11] A.S. Chowdhury, S.M. Bhandarkar, R.W. Robinson, J.C. Yu, Virtual craniofacial reconstruction from computed tomography image sequences exhibiting multi ple fractures, in: Proceedings of 2006 International Conference on Image Pro cessing, IEEE, 2006, pp. 1173–1176, doi: 10.1109/ICIP.2006.312766 .

[12] A.S. Chowdhury, S.M. Bhandarkar, R.W. Robinson, J.C. Yu, Virtual multifracture craniofacial reconstruction using computer vision and graph matching., Com put. Med. Imaging Graph. 33 (5) (2009) 333–342, doi: 10.1016/j.compmedimag. 20 09.01.0 06 .

[13] A. Kronman, L. Joskowicz, Automatic bone fracture reduction by fracture con tact surface identification and registration, in: Proceedings of the 2013 IEEE 10th International Symposium on Biomedical Imaging, IEEE, 2013, pp. 246– 249, doi: 10.1109/ISBI.2013.6556458 .

[14] J. Buschbaum, R. Fremd, T. Pohlemann, A. Kristen, Computerassisted fracture reduction: a new approach for repositioning femoral fractures and planning reduction paths, Int. J. Comput. Assist. Radiol. Surg. 10 (2) (2015) 149–159, doi: 10.1007/s11548 014 1011 2 .

[15] M. Cimerman, A. Kristan, Preoperative planning in pelvic and acetabular surgery: the value of advanced computerised planning modules, Injury 38 (4) (2007) 4 42–4 49, doi: 10.1016/j.injury.2007.01.033 .

[16] A.F. Wilson, P.B. Musgrove, K.A. Buckley, G. Pearce, J. Geoghe gan, D.R. Infirmary, Developing an application to provide interac tive three dimensional visualisation of bone fractures (2009) 10–13. 10.2312/LocalChapterEvents/TPCG/TPCG09/145148 [17] F. Paulano , J. Jiménez , R. Pulido , An application to interact with 3d models reconstructed from medical images, in: Proceedings of the 9th International Conference on Computer Vision Theory and Applications (VISAPP 2014), 2014, pp. 224–229 .

[18] S. Maubleu , C. Marecaux , M. Chabanas , Computeraided planning for zygomatic bone reconstruction in maxillofacial traumatology, in: Proceedings of the 2nd International Conference Surgetica, 2006, pp. 405–411 .

[19] M.H. Moghari, P. Abolmaesumi, Global registration of multiple bone fragments using statistical atlas models: feasibility experiments., in: Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biol ogy Society, vol. 2008, 2008, pp. 5374–5377, doi: 10.1109/IEMBS.2008.4650429 .

[20] R.H. Gong, J. Stewart, P. Abolmaesumi, Reduction of multifragment fractures of the distal radius using atlasbased 2D/3D registration, in: J.P.W. Pluim, B.M. Dawant (Eds.), SPIE Medical Imaging, 2009, pp. 726137–726137–9, doi: 10. 1117/12.811638 .

[21] T. Albrecht, T. Vetter, Automatic Fracture Reduction, in: J.A. Levine, R.R. Paulsen, Y. Zhang (Eds.), Mesh Processing in Medical Image Analysis 2012, Lecture Notes in Computer Science, vol. 7599, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 22–29, doi: 10.1007/978 3 642 33463 4 _ 3 .

[22] T.P. Thomas, D.D. Anderson, A.R. Willis, P. Liu, M.C. Frank, J.L. Marsh, T.D. Brown, A computational/experimental platform for investigating three dimensional puzzle solving of comminuted articular fractures., Comput. Methods Biomech. Biomed. Eng. 14 (3) (2011) 263–270, doi: 10.1080/ 10255841003762042 .

[23] Z. Kato, A unifying framework for correspondenceless shape alignment and its medical applications, Commun. Comput. Inform. Sci. 276 (2013) 40–52, doi: 10. 1007/978 3 642 37463 0 _ 4 .

[24] C. Letta, A. Schweizer, P. Fürnstahl, Quantification of contralateral differences of the scaphoid: a comparison of bone geometry in three dimensions, Anatomy Res. Int. 2014 (2014) 5, doi: 10.1155/2014/904275 .

[25] F. Paulano, J. Jiménez, R. Pulido, Fractured bone identification from CT images, fragment separation and fracture zone detection, in: Developments in Medical Image Processing and Computational Vision, Springer International Publishing, 2015, pp. 221–239, doi: 10.1007/978 3 319 13407 9 _ 14 .

[26] R.B. Rusu, S. Cousins, 3D is here: Point Cloud Library (PCL), Proceedings of 2011 IEEE International Conference on Robotics and Automation, IEEE, 2011, doi: 10. 1109/ICRA.2011.5980567 .

[27] C.C. Ho, F.C. Wu, B.Y. Chen, Y.Y. Chuang, M. Ouhyoung, Cubical march ing squares: adaptive feature preserving surface extraction from volume data, Comput. Graphics Forum 24 (3) (2005) 537–545, doi: 10.1111/j.14678659.2005. 00879.x .

[28] F.R. Feito, J.C. Torres, Inclusion test for general polyhedra, Comput. Graphics 21 (1) (1997) 23–30, doi: 10.1016/S0 0978493(96)0 0 0672 .

[29] P.J. Schneider , D. Eberly , Geometric Tools for Computer Graphics, Elsevier Sci ence Inc., New York, NY, USA, 2002 .

[30] P.J. Besl, D. McKay Neil, A method for registration of 3D shapes, IEEE Trans. Pattern Anal. Mach. Intell. 14 (2) (1992) 239–256, doi: 10.1109/34.121791 .

[31] O. Ron, L. Joskowicz, C. Milgrom, A. Simkin, Computerbased periaxial ro tation measurement for aligning fractured femur fragments from CT: a feasibility study., Comput. Aided Surg. 7 (6) (2002) 332–341, doi: 10.3109/ 10929080209146522 .

[32] S.M. Bhandarkar, A.S. Chowdhury, Y. Tang, J.C. Yu, E.W. Tollner, Computer vision guided virtual craniofacial reconstruction, Comput. Med. Imaging Graph. 31 (6) (2007) 418–427, doi: 10.1016/j.compmedimag.2007.03.003 .