A System for 3D Reconstruction Of Comminuted Tibial Plafond Bone Fractures
CComputerized Medical Imaging and Graphics (2021)
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics
A System for 3D Reconstruction Of Comminuted Tibial Plafond Bone Fractures
Pengcheng Liu a , Nathan Hewitt a , Waseem Shadid a,b , Andrew Willis a, ∗ a University of North Carolina at Charlotte, 9201 University City Blvd., Charlotte, NC 28223-0001, USA b LEAD Technologies, Inc., 1927 South Tryon Street, Suite 200, Charlotte, NC 28203, USA
A R T I C L E I N F O
Article history : Keywords : 3D puzzle solving, computa-tional 3D reconstruction, bone fracturesegmentation, tibial plafond fracture
A B S T R A C THigh energy impacts at joint locations often generate highly fragmented, or commin-uted, bone fractures. Current approaches for treatment require physicians to decide howto classify the fracture within a hierarchy fracture severity categories. Each categorythen provides a best-practice treatment scenario to obtain the best possible prognosis forthe patient. This article identifies shortcomings associated with qualitative-only eval-uation of fracture severity and provides new quantitative metrics that serve to addressthese shortcomings. We propose a system to semi-automatically extract quantitativemetrics that are major indicators of fracture severity. These include: (i) fracture surfacearea, i.e., how much surface area was generated when the bone broke apart, and (ii)dispersion, i.e., how far the fragments have rotated and translated from their originalanatomic positions. This article describes new computational tools to extract these met-rics by computationally reconstructing 3D bone anatomy from CT images with a focuson tibial plafond fracture cases where di ffi cult qualitative fracture severity cases aremore prevalent. Reconstruction is accomplished within a single system that integratesseveral novel algorithms that identify, extract and piece-together fractured fragmentsin a virtual environment. Doing so provides objective quantitative measures for thesefracture severity indicators. The availability of such measures provides new tools forfracture severity assessment which may lead to improved fracture treatment. This paperdescribes the system, the underlying algorithms and the metrics of the reconstructionresults by quantitatively analyzing six clinical tibial plafond fracture cases. ©
1. Introduction
Accurate reconstruction of a patient’s original bone anatomyis the desired outcome for surgical treatment of a bone frac-ture. Treatment goals include achieving expeditious recon-struction and avoiding Post-Traumatic OsteoArthritis (PTOA).When there is involvement of an articulating joint such as thehip, knee, or ankle, accurate reconstruction of the bone joint ∗ Corresponding author e-mail: [email protected] (Pengcheng Liu), [email protected] (Nathan Hewitt), [email protected] (Waseem Shadid), [email protected] (Andrew Willis) surface is critical to avoid PTOA. This task can be quite chal-lenging when dealing with highly comminuted fractures. Thisis due to the fact that often dozens of individual fragments areinvolved and they are sometimes displaced significantly fromtheir original anatomic position and scattered in a complex ge-ometric pattern.This article describes a system for 3D medical image andsurface analysis that enables vital new orthopaedic research toimprove treatment for traumatic limb bone fractures. An em-phasis is placed on the tibial plafond fracture as it is both di ffi -cult to treat and can exhibit wide disparity in subjective fractureseverity evaluation. This article details a system constructedas a novel integration of algorithms that jointly enable virtual a r X i v : . [ c s . C V ] F e b / Computerized Medical Imaging and Graphics (2021)
3D reconstruction of a comminuted bone fracture from a 3DCT image of the fractured limb. While the reconstruction re-sult could be used for pre-operative planning, this article de-scribes the application of this system for the purpose of extract-ing quantitative data for fracture severity classification.Accurately classifying the severity of highly comminutedbone fractures can be challenging for orthopedic physicians andsurgeons. Research in (Marsh et al., 2002; Beardsley et al.,2005; Anderson et al., 2008b) states that accurate determina-tion of the initial fracture severity as a ff orded by fracture sever-ity classification is the single most important prognostic deter-minant of long-term joint health subsequent to trauma. Dueto the importance of fracture severity classification, many re-searchers (Charalambous et al., 2007; Brumback and Jones,1994; Thomas et al., 2007; Anderson et al., 2008a; Beardsleyet al., 2005) have investigated this problem where their com-mon goal is to define methods capable of predicting fractureseverity from quantitative measurements derived from medicalimage data. Yet, no approach to date provides a holistic solutionto this di ffi cult problem as an integrated system.Our system accomplishes this task as a sequence of threesteps:1. Fragment surfaces are extracted from CT images,2. Each fragment surface is further decomposed into anatom-ically meaningful sub-regions,3. Fragments are pieced back together in a virtual space witha puzzle-solving algorithm.Quantitative data is extracted for fracture severity assessment isrecovered as a by-product of the puzzle-solving solution. Thesystem enables extraction of previously unavailable quantita-tive information from fractures cases in terms of the bone frag-ment data to assist physicians in determining the clinical frac-ture severity classification of comminuted bone fractures.Results focus on application of this technology to tibial pla-fond fractures. These fractures typically occur as a result ofhigh-energy trauma such as a ballistic penetration, vehicularaccident, or falls from a height. This article focuses on tibialplafond fracture cases for the following reasons:1. The complex characteristics of this kind of fracture canoften create di ffi culties for physicians in making accurateand reliable fracture severity assessments,2. Tibia fractures often involve the ankle joint which is typi-cally di ffi cult to treat,3. The quality of reconstruction is a critical factor for goodprognosis,4. PTOA is directly related to the accuracy of reconstruction.Hence, this sub-domain of fractures represents one of the mostpromising applications of the research described herein.
2. Related Work
Computational puzzle solving in 3D seeks to use computeralgorithms to facilitate reconstructing 3D broken objects fromthe geometry of their fragments. In general, puzzle-solving ap-proaches fall into three categories: (1) boundary matching, i.e., algorithms that match together fragments by comparing theirboundaries, (2) template matching, i.e., algorithms that matchfragments into an a-priori known template that is used as a ref-erence shape for the broken fragments and (3) manual recon-struction approaches. Approaches from the first two categoriesrequire algorithms for curve and surface matching and surfacealignment. For boundary matching, boundary curves from thebroken fragments are matched to piece together the fragments.For template matching, surfaces on the fragment are matchedto the corresponding surfaces on the template so that fragmentscan be aligned into the template to accomplish reconstruction.Boundary matching approaches for puzzle solving starts withseminal work (Freeman H, 1964; Wolfson et al., 1988) that de-veloped algorithms to piece together 2D jigsaw puzzles similarto that depicted in figure 2. However, solving 3D bone frac-ture “puzzles” is a significantly more challenging problem. Jig-saw pieces share similar sizes and have distinctly identifiableshapes that greatly restrict the collection of potential matingsurfaces / curves on corresponding pieces. Further, for jigsawpuzzles one can assume that all of the jigsaw pieces are com-plete and intact, and no pieces are missing. Puzzle-solving3D fractures includes low-resolution sensed data (segmentedCT), deformable fragments (bone tissue), potentially nonde-script fracture surfaces (similar to oblique planes) and due toa high-energy impact there may be small, missing or unusablefracture pieces.Fracture reconstruction from boundary data matches bound-aries and surfaces generated due to the fracture event. Somework on bone fragment reconstruction uses simulated data bysimulating fractures via computer-generated fractures or by us-ing a drop-tower to break bone surrogate material. Examples in-clude (Kronman and Joskowicz, 2013) which simulates a frac-ture by slicing the intact CT scan of a healthy bone. This allowsquick experimentation on a variety of test cases. However, thesesimulations do not capture the complex physiological or biome-chanical phenomena that give rise to these fractures which tendto create fractures that may never been seen in clinical practice.Further, the approach analyzes the geometry of the fracture sur-face generated by separating bone fragments and assumes thesesurfaces will accurately fit together using some fracture surfacematching metric when in reality bone fragment fracture surfacesmay not geometrically match well. In (Willis et al., 2007) tibialplafond fractures are simulated using a drop tower to simulatehigh axial load to the ankle joint to fracture surrogate bone con-structed from high-density polyether urethane foam. Result-ing fragments were subsequently scanned in a laboratory set-ting an fracture surfaces were matched to puzzle-solve the frac-ture. Here, the measured data and extracted surfaces exhibitlittle noise which makes reconstruction much easier. Work in(Fürnstahl et al., 2012) proposes an approach for reconstructingproximal humerous joint by alignment of fracture surfaces. Aconcern here is the geometric accuracy of extracted fracture sur-face segmentation; especially in low-contrast trebecular bonetissue regions. We opt to align outer bone surfaces only avoid-ing the potentially inaccurate measurements of the fracture sur-face regions.A second category of approaches for bone fragment recon- Computerized Medical Imaging and Graphics (2021) 3Figure 1: A brief overview of a system proposed in this paper. This system takes as input a 3D CT image of the fractured limb and a 3D CT image of the undamaged(intact) limb, and provides as output a virtual reconstruction of bone fragments which estimates the anatomy of the patient’s original bone. Each block denotes astep and the associated algorithm name is shown inside the block (in blue). Images show how the data has been changed in each step of the 5-step reconstructionprocess. (a) (b)Figure 2: The di ff erence between reassembling (a) commercial jigsaw puzzlesand reconstructing (b) broken artifacts. Note here that the latter is made muchmore di ffi cult as corners are not easily identifiable and may not indicate thebeginning or end of a curve that will uniquely match some other fragment.Worse, any portion of a boundary curve may match to any other fragment andthe curve itself may match equally well with numerous similar boundaries fromother fragments. (Used with permission from Willis and Cooper (2008). ) struction focuses on aligning fragments into a template of theintact bone. Alignment is often achieved through applicationsof the Iterative Closest Point (ICP) algorithm. This has thedistinct disadvantage that, especially in cases where fragmentshave undergone large displacements or rotations, the ICP algo-rithm can converge on local minima, giving inaccurate results.In (Chowdhury et al., 2009), work is done on reconstructing fa-cial fractures. As the geometry here is more unique than thetibia, the assumption can be made that each bone surface willhave only one match. Because of the di ff erences between com-minuted tibia fractures and craniofacial fractures, the e ff ective-ness of ICP in this context translates poorly to our application.The work done in (Albrecht and Vetter, 2012) demonstratesthat, with modification, the ICP algorithm can improve uponcurrent reconstruction standards. While a deliberate e ff ort ismade to lessen the likelihood of an erroneous convergence ofthe algorithm, it is still possible and is thus an undesirable algo-rithm. Other work in (Okada et al., 2009) reconstructs frag-ments of the proximal femur for fracture cases involving upto 4 fragments by matching fragment outer surfaces and theirboundaries. While multiple fragment solutions are described, there are examples of a multi-fragment solution and no compu-tational framework to cope with the large number of potentialmatches between fragments, fragment boundaries and their un-known potential locations. (Moghari and Abolmaesumi, 2008)defines an Euclidean invariant coordinate system within whichto fit algebraic, i.e., implicit polynomial, surface representationsto bone fragments surfaces and uses the coe ffi cients of the poly-nomial models to identify matches between fragments and theircorresponding atlas locations.Several bone fragment reconstruction methods have beendocumented that are dependent on manual input. In (Scheuer-ing et al., 2001), a system is proposed that simulates volumet-ric collision detection in a virtual 3D environment and an op-timization process for repositioning the bone fragments. Inthis application, users manually place the fragments close to-gether, then refine the alignment via computer optimization.This method is made more interactive by the work in (Harderset al., 2007), which takes a similar approach of user-guided re-construction with collision detection while adding haptic feed-back to the system. This more advanced human-computer in-teraction increases the intuitiveness of the system for users.Work in (Cimerman and Kristan, 2007; Kovler et al., 2015) pro-poses user-guided reconstruction with collision detection whileadding surgery simulation. These reconstruction approaches allfocus on manual, i.e., interactive alignment, whose applicationsinclude pre-operative planning, virtual surgical procedures orsurgical training. In contrast, this work targets objective extrac-tion of fracture metrics for severity analysis. As such, we seekto minimize interaction which can bias computed results.This work represents a significant advancement in fracturereconstruction state-of-the-art in several important ways:1. it demonstrates reconstruction solutions for clinical datawhich include the most comminuted, i.e., complex, frac-ture cases in the literature to date (10 fragments).2. it incorporates computational acceleration methods re-quired by complex cases which significantly reduce com-putational cost. / Computerized Medical Imaging and Graphics (2021)
3. it exhibits smaller average reconstruction alignment errorsthan observed in competing solutions (Okada et al., 2009;Moghari and Abolmaesumi, 2008; Fürnstahl et al., 2012).Further, the end-purpose of the system as a means to facilitatefracture severity assessment is novel relative to other implemen-tations which focus on surgical planning.
3. Methods
Figure 1 depicts the proposed approach for fracture recon-struction as a sequence of five steps. Implementation requiresfive application-specific algorithms that process 3D image andsurface data. These five algorithms are:1. 3D CT image segmentation (see §3.3)2. 3D surface partitioning (see §3.4)3. Appearance-based 3D surface classification (see §3.5)4. 3D puzzle-solving (see §3.6)5. Quantitative fracture analysis (see §3.7)A detailed list of processing steps, including the intact templateprocessing is described in Algorithm 1. These steps also specifythe sequence of algorithms through which the data flows as eachalgorithm requires results computed from the previous one.
We provide a computational complexity analysis of our en-tire system by assessing the complexity of each of the stepsrequired to compute the reconstruction, i.e., steps 1-4 above.The segmentation algorithm of step (1) has a computationalcomplexity O ( M + N log( N )) Felkel et al. (2001) where N de-notes the total number of voxels in the medical image and M denotes the number of edges between voxels as defined by theadjacency scheme. The algorithm we use, (Shadid and Willis,2018, 2013), implements a 6-adjacency scheme which M (cid:39) N and the resulting algorithm have complexity O ( M + N log( N )).The surface partitioning algorithm of step (2) uses the graphof the surface mesh to partition a closed surface into surfacepatches. Surface partitions are defined using a minimum span-ning tree of the surfaces which has computational complexity O ( E log( V )) for a mesh surface having E edges and V vertices.The classification algorithm of step (3) requires O ( N ) calcula-tions to complete. The 3D puzzle solving algorithm of step (4)consists of two parts: (i) a surface correspondence search and(ii) candidate match evaluation. The surface correspondencesearch problem is solved via brute-force pairwise evaluation of N surface correspondences. Each correspondence is assigneda similarity score by matching spin image surface features, analgorithm that has linear computational complexity in the spinimage descriptor dimension, i.e., O ( N overlap ), where N overlap de-notes the number of overlapping pixels in a spin image pairmatch (see § 3.6.3). A small subset of surface matches havinghigh similarity scores are processed using the ICP algorithm tocompute the final fracture reconstruction which has computa-tional complexity O ( kV o ) where V o denotes the number of sur-face vertices and k denotes the number of iterations required forICP to converge. We mention k explicitly to point out that our method for coarsely aligning surfaces via spin images can sig-nificantly reduce k by starting this nonlinear optimization closerto the desired minimum of the error functional. Our system can reconstruct fracture cases completely auto-matically but include interfaces to specify algorithm parame-ters. For some steps, i.e., steps (4) and (6) of from Algorithm1 an interface has been included to allow manual interaction tosteer the system to a better solution. Step (4) from Algorithm1 requires access to a CT-machine specific training databaseto classify fragment surfaces. If this data is not available, aninterface is provided that allows the user to manually classifyfragments surfaces with may take several minutes to specify in-teractively. Step (6) from Algorithm 1 also has an interactiveinterface which allows users to manipulate fragments in the fi-nal solution. This interaction is rarely used and typically is ap-plied to re-initialize fine alignment, e.g., one may slightly rotatea small fragment and initiate new a global puzzle alignment so-lution to find a better solution.
Reconstruction of a 3D solid from its fragments is a geo-metric problem where the geometry of the fragments must beknown in order to piece them back together. For this reasongeometric models for each bone fragment in the 3D CT imagemust be computed. This is the goal of the CT image segmen-tation process. The segmentation algorithm used assigns eachpixel in the CT image I ( x , y , z ) to a label, l . The goal of the sur-face segmentation algorithm is to estimate the correct label forevery pixel and we denote this operation is mathematically with F { I } . Fragment surfaces may be extracted from the labeled im-age by estimating the locations where the image changes label,or equivalently, solving for the locations where the label im-age changes from fragment label, l = k , to a di ff erent labelvalue, i.e., S k = (cid:26) F (cid:110) I f (cid:111) | (cid:13)(cid:13)(cid:13)(cid:13) ∇F (cid:110) I f (cid:111)(cid:13)(cid:13)(cid:13)(cid:13) (cid:44) , F (cid:110) I f (cid:111) = k (cid:27) . Geo-metric segmentation of bone fragments is challenging due tothe similarity in cancellous and background CT tissue inten-sities making standard (van Eijnatten et al., 2018) or manual(Paulano et al., 2014) approaches impractical. This work adoptsthe approach described in (Shadid and Willis, 2018, 2013) tosegment the image which improves upon prior watershed al-gorithm implementations (F.Meyer, 1994); especially for bonefragment segmentation. Bone fragment surfaces are extractedfrom the segmented image using the marching cubes algorithm(Lorensen and Cline, 1987). There are two major challenges tothis segmentation problem:1. Bone fragment boundaries are especially di ffi cult to de-marcate when fragments are abutting other fragments.2. Intensities for some bone tissues, i.e., cancellous tissue,are the same as that for other soft tissues in the sameanatomic region making them di ffi cult to discriminate.The bone fragment segmentation algorithm takes as input a CTimage of the limb and two user parameters: (1) the maximum Computerized Medical Imaging and Graphics (2021) 5
Algorithm 1
Reduction of bone fragment k in 7 algorithmic steps. Processing requires (5) algorithms: (1) an CT image bonesegmentation function, F { I } , (2) a fragment surface segmentation function, F {S} , (3) a fragment surface classifier,
C {S} , (4) asurface feature descriptor extractor,
E {S} and (5) a puzzle solver, P (cid:110) D k , o , D T (cid:111) . Lines 3-7 are repeated until all fracture fragmentsare aligned. Steps ending with an asterisk ∗ above have facilities to assist the solution via manual interaction (see § 3.2 for details). Input:
Two 3D CT images: (1) a fractured limb CT, I f ( x , y , z ), and (2) an intact (contralateral) limb CT , I i ( x , y , z ). Output:
The reduced fracture poses, T k , and fracture severity analysis data, A k .1: S T = F { I i } : The template bone surface S T is extracted from the intact CT image I i . O ( N )2: D Tj = E (cid:110) S T (cid:111) : Surface descriptors, D Tj , are extracted from the template bone surface, S T .3: S k = F (cid:110) I f (cid:111) : The k th bone fragment surface S k is extracted from the fracture CT image I f .4: S k , o = C (cid:110) F (cid:110) S k (cid:111)(cid:111) : Extract the k th fragment outer / cortical surface, S k , o , by partitioning and classifying the fragment surface, S k . ∗ D k , o = E (cid:110) S ko (cid:111) : Extract a set of surface descriptors, D k , o , from the k th fragment outer / cortical surface, S k , o .6: T k = P (cid:110) D k , o , D T (cid:111) : Estimate the k th fragment’s pose, T k , from matched template / fragment feature pairs, (cid:110) D k , oi , D Tj (cid:111) . ∗ A k = L { T k } : Fracture severity analysis data, A k , is extracted from the reduced fracture solution and surface data, (cid:110) T k , S k , o (cid:111) . (a) (b) (c) (d)Figure 3: (a) shows slices from a 3D CT image of an intact ankle joint. (b)shows the segmented tibia, fibula, and talus bone surfaces. (c) shows slices froma 3D CT image of a fractured ankle. (d) shows the segmented bone fragmentssurfaces. The segmentation from the 3D CT image was performed using theproposed segmentation algorithm. size of a bone region, and (2) a segmentation sensitivity param-eter. The output of this algorithm is a labeled CT image whereeach unique label corresponds to a unique bone fragment.The segmentation algorithm proceeds by incrementally clas-sifying pixels from image regions of high-likelihood (corticaltissue) to regions of low-likelihood. To do so, image pixelsare initially classified into three sets: (1) non-bone, (2) corti-cal bone, and (3) non-cortical bone. Connected cortical bonepixel regions are assigned a fragment label and a customizedwatershed algorithm (Shadid and Willis, 2013), referred to asthe Probabilistic Watershed Transform (PWT), merges adjacentpixels into each fragments region. The merge procedure calcu-lates the conditional probability of each candidate pixel giventhe current segmentation and includes models to specifically ad-dress the realization of layered, i.e., lamellar bone tissue struc-tures in CT images, and common fracture phenomenon suchas fragment splintering which generates long, thin protrusionsof bone tissue that is di ffi cult to otherwise segment. Surfacesare extracted from the segmented image data with the marchingcubes algorithm (Hansen and Johnson, 2005). Figures (3(b,d))show 3D surfaces segmented from CT images of ankle joints foran intact (b) and a fractured (d) ankle joint using the proposedalgorithm. Puzzle-solving requires shapes to be decomposed into shapeparts. Similarly, for fracture reconstruction, each segmentedbone fragment model must be partitioned into a collection ofsurface patches to allow the outer surfaces of each fragment sur-face to be geometrically matched to the template surface, S T . Asurface segmentation algorithm divides fragment surfaces into3 anatomically distinct groups which greatly reduces the di ffi -culty of the search problem by limiting the number of plausiblecorrect surface matches and thereby improving the reconstruc-tion system performance. The surface mesh partitioning algo-rithm divides the k th bone fragment surface, S k , into a collec-tion of surface patches, { s , s · · · } , such that S k = ∪{ s , s · · · } .Each of the generated surface patches are intended to consist ofsurface points from only one of three anatomical categories. Todo so, the system applies a “ridge walking” algorithm (Willisand Zhou, 2010) which solves this problem geometrically us-ing the fact that the anatomical categories of interest can bediscriminated well by dividing the 3D bone fragment surfacesalong 3D contours that traverse high-curvature ridges. A surface classification algorithm must identify the outer,i.e., cortical, surface(s) of each bone fragment from the col-lection surface patches created by the surface partitioning al-gorithm. This is a classification problem that is solved usingby a classification algorithm,
C { s , s · · · } that extracts the pe-riosteal, i.e., outer, surface, S k , o , from the set of partitioned sur-faces, { s , s · · · } from the k th fragment for template-based re-construction. Our novel classification approach uses CT bonetissue intensity variations in the vicinity the fragment surface todetect the unknown anatomic label of the surface. This is madepossible by the fact that cortical, trebecular, cancellous and sub-chondral bone tissues have distinct intensity ranges and relativethicknesses in the anatomic regions of interest. Unfortunately,intensity variations by machine, patient gender and age (Lookeret al., 2009; Felson et al., 1993) require the user to specify train-ing data from within the CT image for reliable output. Figure4 shows image data from the three anatomic bone regions thatmust be marked interactively by a user: (1) the diaphysis, made / Computerized Medical Imaging and Graphics (2021)Figure 4: CT appearance of the distal tibia anatomy. Relative to the tibia’s outersurface, the CT intensities vary with characteristic profiles along the inwardsurface normal for each anatomic region.Figure 5: Two views of the surface mesh classification results for two frag-ments. Red denotes periosteal surface, green denotes fracture surface, and bluedenotes articular surface. of solid dense cortical bone, (2) the metaphysis, made of a cor-tical shell and an interior, porous, cancellous bone, and (3) theepiphysis, made of dense subchondral bone. From this data, theautomated classification algorithm (Liu, 2012) assigns the parti-tioned fragment surface patches to one of three semantic labels:(1) “fracture surfaces” (surfaces generated when the bone brokeapart), (2) “periosteal surfaces” (surfaces that were part of theouter bone surface), and (3) “articular surfaces” (surfaces thatfacilitate the articulation of joints). Classification of these sur-face patches is useful for both bone reconstruction and severityanalysis. For reconstruction, the periosteal surfaces for geomet-ric alignment for fragments template-based reconstruction (theapproach used in this paper), (2) fracture surfaces allow calcu-lation of surface area produced by a fracture event, which isa valuable objective severity metric, and (3) articular surfacesare particularly medically important in reconstruction. Fromthis data a classifier is constructed and applied to automaticallyclassify the extracted surface patches to these three anatomiclabels. A final step merges adjacent periosteal surface patcheshaving the same anatomic label generating the periosteal, i.e.,outer, surface, S k , o . Figure 5 shows classification results for twobone fragments using this approach. The 3D puzzle-solving algorithm takes as input the k th frag-ment outer surface patch, S k , o , and estimates the transforma-tion, T k , that rotates and translates this fragment to it’s originalanatomic position within the intact template. A complete virtualreconstruction of the unbroken bone from the bone fragments isachieved by transforming the geometry of each fragment usingit’s associated transformation. The intact contra-lateral bone, asrepresented in the intact CT image, is taken as a reasonable ap-proximation of the unbroken bone after mirroring the geometryacross the human bilateral plane of symmetry. Puzzle solving algorithms solve a di ffi cult computationalsearch problem whose goal is to estimate the unknown geomet-ric correspondences between puzzle pieces. This is typicallyaccomplished in two steps: (1) hypothesize geometric matchesand (2) test hypothesized matches by evaluating a test statistic;typically the minimized pairwise surface alignment error usingthe ICP algorithm (Besl, 1992; Chen and Medioni, 1991; Low,2004; Zhang, 1994)). Each hypothesis has the form: “Fragmentsurface, S k , o , corresponds template surface, S T , at matchingsurface points { p k , p T } respectively.” For the k th puzzle frag-ments the puzzle-solved fragment pose, T k , is taken as the pair-wise hypothesis which, after alignment, has the smallest align-ment error.Use of the ICP algorithm to puzzle solve the fracture is bothan unreliable and computationally prohibitive option. Reliabil-ity is an issue due to the tendency of the ICP algorithm to con-verge to the closest local minima. Hence, unless the fragmentsare close to their final anatomic pose in the puzzle solution, itis unlikely their final alignment will be correct. Computationalcost is an issue since many fragment surfaces include thousandsof samples and finding corresponding surface points is a O ( N )search problem and evaluation of each pairwise correspondencerequires a iterative ICP solution which has computational com-plexity O ( N ). Hence the total computational cost of ICP-basedreconstruction for one fragment is O ( N ) where N is the num-ber of surface points. For m fragments the complexity increasesto O ( mN ) ( m (cid:28) N ) (Arthur and Vassilvitskii, 2006). For ourexperiments, N ∼ which imposes computational costshigh enough to make brute force solutions computationally pro-hibitive for even very small puzzles. An approach for bone fragment puzzle solving is proposed tomake the puzzle solving algorithm computational cost tractable.Descriptor-based puzzle solving approaches extracts descrip-tors from surfaces and uses these descriptors to hypothesizematches and geometrically align the pieces and use the averagealignment error as a statistic to evaluate the likelihood that thehypothesized match is correct. The approach to puzzle-solvethe k th fragment, S k , o , with respect to intact template, S T , con-sists of five steps:1. Descriptor Extraction : This step applies a surface de-scriptor extractor algorithm,
E {S} , to puzzle surfaces, S ,to produce a set of Euclidean-invariant surface descrip-tors. Let D T = E (cid:110) S T (cid:111) denote descriptors extracted fromthe template surface and D k , o = E (cid:110) S T (cid:111) denote descriptorsfrom the k th fragment surface. (see 3.6.3)2. Generate Matched Features:
This step takes in as in-put the computed descriptors from (1) and outputs a list ofmatched descriptor pairs: L initial . (see 3.6.4)3. Remove Incorrect Matches:
This step takes L initial as in-put and removes suspected false matches from the initiallist and outputs a list of candidate descriptor matches in anew list: L candidate . (see 3.6.5) Computerized Medical Imaging and Graphics (2021) 7 Test Candidate Matches:
This step takes in the previ-ously generated list L candidate and outputs a 3D transfor-mation matrix, T k , for each match which aligns the bonefragment to the surface of the intact template.(see 3.6.6)5. Select The Best Matches:
This step determines which ofthe alignments from the previous step provide the best re-sult and uses these to provide the final puzzle-solved solu-tion. (see 3.6.7)Applications in surgical planning and image guided surgery re-quire the puzzle solution to be geometrically registered to thecoordinate system of the fractured CT image, I f . To this end,an initial alignment is performed between intact template sur-face, S T , and the anatomically un-perturbed fragment of thefractured limb which we refer to as the “base fragment.” Fortibial plafond fractures the “base fragment,” is typically the up-permost (proximal) bone fragment in the fracture. This initialmovement of the template serves to bring the template surfaceinto the coordinate system and vicinity of the bone fragmentsurfaces improving convergence rates of iterative geometricalignment optimizations. Figure 6 shows a graphical overviewof this method for puzzle-solving a clinical bone fracture case. While there are many candidate surface descriptors that couldbe used to compactly characterize local surface patch geometry,this paper adopts the spin image representation (Johnson andHebert, 1997, 1996) as its surface descriptor. The rationale forthis choice is based on the low computational cost this repre-sentation a ff ords when performing surface matching. Specifi-cally, spin images converts a nonlinear 3D surface-pair align-ment problem into a low-dimensional 2D (spin) image-paircomparison problem (see figure 7d). Spin images of two iden-tical surface patches in arbitrary pose generate produce iden-tical 2D spin images. As mentioned in § 3.1, use of spin im-ages significantly reduces the complexity of comparing surfaceshapes by reducing the computational complexity from poly-nomial time O ( kN )(ICP) to linear O ( N overlap ) for spin imagepair having N overlap overlapping pixels (2D-image subtraction).Unfortunately, matching spin images does not solve the com-plete alignment problem and requires estimation of a singleadditional parameter: rotation about the corresponding surfacenormals. Hence, spin images provide a compact, yet imperfect,Euclidean invariant representation of the shape of a bone frag-ment surface. Let L initial denote the list of initial hypothesized matcheswhere each match is specified by an index pair { i , j } indicat-ing that the i th surface descriptor from the template, D Ti , is hy-pothesized to correspond to the j th surface descriptor from the k th fragment’s outer surface D k , oj . Evaluation of the correspon-dence hypotheses is accomplished by calculating the value of asimilarity function C sp ( P , Q ) to detect similar surface matches.To avoid biases in descriptor matching, the proposed functioncombines statistical correlation and sample size into this simi-larity function. Our similarity function uses the normalized linear correla-tion coe ffi cient, R sp ( P , Q ), to measure descriptor similarities asshown in equation ((1)) which assigns a score of 1 for a perfectmatch and scores can range from (-1, 1). Given two spin im-ages P and Q with N bin bins (number of pixels in the image ),denote p i = P ( α, β ) as the i th spin image intensity value fromspin image P , similarly, denote q i = I ( α, β ) as the i th spin im-age intensity value from spin image Q . With this notation, thelinear correlation coe ffi cient, R sp ( P , Q ), is computed as shownin equation (1): R sp ( P , Q ) = N bin (cid:80) i p i q i − (cid:80) i p i (cid:80) i q i (cid:115) N bin (cid:80) i p i − (cid:32)(cid:80) i p i (cid:33) N bin (cid:80) i q i − (cid:32)(cid:80) i q i (cid:33) (1) When R sp is close to 1, the images are similar, when R sp isclose to -1 the images are di ff erent. Unfortunately, this met-ric is biased since each spin image may have di ff erent ( α, β )dimensions. Hence good matches may be found when spin im-ages have similar values over small overlapping ( α, β ) regions.On the other hand, two spin images with a large area of over-lap may be assigned a lower R sp value but exhibit widespreadsimilarities.To address biases R sp ( P , Q ) the similarity function C sp ( P , Q )proposed in (Johnson and Hebert, 1997, 1996) was used asshown in equation (2). This metric is well-suited because itconsiders both the correlation, i.e., match fidelity, and also thenumber of matching or overlapping pixels, N overlap , as part ofthe final matching score. C sp ( P , Q ) = (cid:16) atanh (cid:16) R sp ( P , Q ) (cid:17)(cid:17) − λ (cid:32) N overlap − (cid:33) (2)The resulting similarity function weights the correlation R sp against the variance intrinsic to the correlation coe ffi cient whichincreases as the amount of overlap in the two spin images in-creases. This is accomplished via a change of variables com-monly referred to as "Fisher’s z-transformation" (Fisher, 1915)which uses the hyperbolic arctangent function to transform thedistribution of the correlation coe ffi cient into an approximatelynormal distribution having constant variance N overlap − . This al-lows the similarity metric (2) to be formed that strikes a com-promise between good correlations (term 1) and su ffi cient ev-idence, i.e., measurements to trust the computed correlation(term 2) to better detect reliable matches. The parameter λ isused to control the relative weight of the expected value of thecorrelation coe ffi cient and the variance of this statistic to pro-duce a final similarity score. (Johnson and Hebert, 1997, 1996)mention that λ controls the point at which the overlap betweenspin images dominates the value of the similarity metric for twospin images. When the overlap is much larger than λ , the sec-ond term in equation (2) becomes negligible. In contrast, whenthe overlap is much less than λ , the second term dominates thesimilarity measure. Therefore, λ should be the expected overlapbetween spin images. In this paper, λ is automatically computedfor each surface match by computing average number of non-black pixels for all the spin images generated from that frag-ment and setting λ to half of the average value. / Computerized Medical Imaging and Graphics (2021)Figure 6: Puzzle solving process: initialization, matching and alignment.(a) (b) (c) (d)Figure 7: (a) Oriented surface point p with its normal n and the tangent plane p , one of its neighbor point x on the fragment surface and ( α , β ) are computedvalues for point x respect to point p ; (b) shows point p on the bone surface andits neighbor points; (c) Plotted neighbor points on the ( α, β ) grid; (d) Computedspin image for out of (c) Johnson and Hebert (1997). Due to the noise from image data and errors from segmen-tation and classification algorithms, the initial list of matches L initial often contains many false hypothesized surface point cor-respondences. Since careful analysis of these matches is com-putationally expensive, we propose an alternative approach todetect incorrect hypotheses that significantly improves perfor-mance. To do so we consider multi-hypothesis consistency con-straints to e ffi ciently eliminate implausible solutions.Conceptually, the idea here is that each fragment has a singleunique Euclidean transformation that brings the fragment to it’soriginal anatomic location, T k . Hence, if multiple hypothesesfrom the same fragment are made, those that are correct musthave nearly the same value for T k . Towards this end, evalu-ated hypotheses for the k th can be used to cross-validate otherhypotheses from the same fragment. However, rather than per-forming the computationally expensive alignment required tocompute T k , we use similarity as provided by spin image sur-face descriptors to reject large numbers of incorrect hypotheses.Let (cid:110) D Ti , D k , oj (cid:111) and (cid:110) D Tl , D k , om (cid:111) denote descriptor pairs { i , j } and { l , m } from the list L initial computed from surface points are (cid:110) p Ti , p k , oj (cid:111) and (cid:110) p Tl , p k , om (cid:111) . If both hypotheses are correct, then thesurface point pairs should be geometrically consistent, i.e., thedistance between p k , oj and p k , om should be equal to the distancebetween p Ti and p Tl . This constraint can be directly validatedusing spin image coordinates in a geometric consistency test.If the two matches satisfy equation ((3)), we consider them tobe geometrically consistent matches which means both of them may be true matches. (cid:12)(cid:12)(cid:12)(cid:12) D Ti (cid:16) p k , oj (cid:17) − D k , oj (cid:16) p Ti (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12)(cid:12) D Tl (cid:16) p k , om (cid:17) − D k , om (cid:16) p Tl (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) < D gc (3)Where D gc = γ intact where γ intact is the resolution of theintact template, i.e., the average edge length of the edges fromthe intact template surface mesh.When the initial match list is constructed, list elements aresorted by decreasing similarity score, C sp . Inconsistent ele-ments are removed by splitting the list L initial into two partsat its midpoint. The first list contains matches having highersimilarity scores and the second contains matches that havinglower similarity scores. We consider one match from eachof the two lists respectively form pairwise correspondence hy-potheses. The geometric consistency condition is evaluated forthe pair of correspondences. If they both satisfy the consis-tency condition, we keep both matches. Otherwise, we keep thematch that has a higher similarity score and discard the othermatch. After evaluating all matches in both lists, the remainingmatches are placed in the sorted list referred to as L candidate . The reduced set of hypotheses in the L candidate set are nowevaluated by calculating their 3D surface alignment error. Asmentioned previously, this is a computationally expensive stepinvolving nonlinear optimization via the ICP algorithm whichis further known to su ff er from erroneous solutions when theinitial transformation, T , for optimization lies far from the truevalue, T k . Our surface descriptor matches serve to significantlyreduce computation here by performing a coarse initial align-ment which provides a good initial guess for 5 of the 6 unknowntransformation parameters, specifically the (3) translation pa-rameters and (2) of the three unknown rotation parameters. Fora given match, the alignment proceeds in two steps:1. Coarse Surface Alignment: For each likely descriptormatch hypothesis (cid:110) D Ti , D k , oj (cid:111) an initial transform, T , iscalculated which translates the fragment surface, S k , o , tomake it’s surface point, p k , oj , correspond with the tem-plate’s surface point p Ti and simultaneously rotates thefragment surface to make the fragment and template sur-face normals also correspond at these points. Computerized Medical Imaging and Graphics (2021) 9Figure 8: An example of alignment using the surface alignment engine. Thered bone is the intact template, the green bone is one of the fracture fragment.The two black points on both surfaces are matched surface points.
2. Refined Surface Alignment: Using the initial transform T , we run the ICP algorithm to compute final 3D align-ment error.In this way hypotheses can be used that may not be exactlycorrect and in cases where the hypothesis is “close,” the ICPalgorithm is still likely to converge to the correct alignment. Figure (8) illustrates how the alignment process uses a hy-pothesized surface correspondence to align a fragment to theintact template. For a given fragment, each hypothesized sur-face correspondence is evaluated in the sequence given by thematch score from §(3.6.5) starting at the highest score. Eachevaluation aligns a fragment to the intact template using thethree alignment steps specified in the three previous section.However, if the algorithm goes through all three steps for ev-ery match, the reconstruction process is very time-consuming.To reduce computation each alignment step includes predefinedconditions to determine whether the system should continue toevaluate the match or discard the match and try a new hypoth-esized match. The local alignment error is used to control theprocess. Because the computational cost increases in each step,the alignment error threshold values are smaller (stricter) foreach step. In this paper, the threshold values are 4 γ intact for stepone, 2 γ intact for step two, and γ intact for step three. Finally, whenone match is accepted by all three steps the output position fromstep three is considered as the final alignment for the fragment. Several software enhancement tools are introduced to helpimprove the quality of the 3D puzzle-solved solutions and tohelp reduce the time necessary to compute these solutions.These enhancement tools include: (1) surface sampling, (2) us-ing occupied regions, (3) Mean Curvature Histogram BiasedSearch, and (4) global fragment alignment optimization. Thefollowing sections (3.6.9, 3.6.10, 3.6.11, 3.6.12) discuss eachtool in detail.
The intact template and the fragment surfaces often consistsof a large number of surface points. Computing spin images forevery point on these surfaces is a time-consuming task. In or-der to improve the speed of the puzzle-solving algorithm, uni-form sub-sampling is applied on both the intact template andthe fragment surfaces. Surface points are randomly selected on the fragment with a constraint that the distance between anytwo sampled surface points are greater than a predefined sam-pling distance (cid:52) s . The larger (cid:52) s value computes fewer sam-pled points on the fragment surface and smaller value com-putes more sampled points on the fragment surface. Spin im-ages are only computed for those sampled surface points onthe intact template and the fragment surfaces. In this paper,the sub-uniform sampling distance for each fragment is set to (cid:52) s = . γ intact , and γ intact is the average edge length of edgeson the intact template. The concept of occupied regions allows for significant per-formance improvements by further reducing the spin imagescomputed on the intact template. Since part of the intact tem-plate surface has already been aligned to the base fragment,these points should be excluded from the matching. We markthese surface points on the intact template surface as “occu-pied”, and spin images are only computed for surface pointsinside the “unoccupied” regions of the intact template. Thissignificantly reduces the search space for the matching process.Occupied points on the intact template are flagged using thefollowing condition: the intact surface point, v tj , is marked as“occupied” if its distance to the closest aligned fragment sur-face is less than 2 γ intact . To reduce the number of spin images computed for boththe intact template and fragment surfaces, an approach calledmean curvature histogram biased search is used. Here we biasthe search to emphasize regions having discriminative shapeby selectively computing descriptors in locations where themean surface curvature is large, i.e., highly curved surface lo-cations. The goal is to focus on calculating surface descriptorsthat are more easily matched and have fewer candidate corre-spondences. Using these locations allows computation to fo-cus on regions that carry more information regarding the un-known puzzle solution. Towards this end, each descriptor isattributed with an estimate of the surface mean curvature at theassociated surface point. We then randomly select points uni-formly from the distribution of observed curvatures which dras-tically reduces the number of spin images computed. Further,when forming pairwise hypotheses we only create hypothesesfor points pairs whose mean surface curvature are similar. Thisapproach greatly reduces the computational cost of the puzzlesolving algorithm and improves performance of the reconstruc-tion without an observed sacrifice in solution accuracy.
Global alignment seeks to simultaneously optimize all un-known transformation parameters of the puzzle-solved frag-ments. This serves to correct minor discrepancies in the po-sition of misaligned fragments by fitting together adjacent frag-ments onto the template surface in the final reconstruction re-sults. This is accomplished by perturbing fragment transforma-tion values away from the current puzzle solution and perform-ing ICP simultaneously using the global alignment algorithm / Computerized Medical Imaging and Graphics (2021) of (Pulli, 1999). This approach seeks to simultaneously equal-ized the surface alignment error equally across all fragmentsand enforce the restriction that fragment-and-template surfacecorrespondences are unique, i.e., each template surface pointcan only correspond with points from a single fragment (as de-scribed in 3.6.10). Since the perturbations are small, this pro-cedure serves to make small adjustments each fragments finalsurface position on the intact template and minimize fragment-to-fragment surface overlap. This process can improve the poseof slightly misaligned fragments in the final puzzle solved so-lution.
The post-reconstruction analysis is the final output of the sys-tem. The analysis tools are integrated into the system and allowusers to analyze the reconstruction result and help users betterunderstand the fracture case. The analysis consists of two com-ponents: (a) analysis of the geometric accuracy of the alignedfragments and (b) analysis of the severity of the bone fracture.The first component provides a table containing alignment in-formation, such as sampled, matched, and unmatched points,and a histogram of alignment error for each fragment. Thesecond component is a fracture severity report which containsquantitative values for several key factors for each fragmentwhich are known to be indicators of fracture severity, such asfracture surface area. Although visual assessment of 3D recon-struction result is valuable to users, these analysis tools providequantitative information which can help users objectively inter-pret the fracture case.Analysis of the geometric accuracy of the reconstruction isaccomplished by measuring the point-to-plane distance fromeach vertex of the fragment surface to the plane of the clos-est triangle on the intact template. These distances are aver-aged to compute alignment error for each fragment. The post-reconstruction analysis report includes a table of alignment in-formation and a table of histograms that show the distributionof alignment errors for each fragment.Our system provides histograms of the alignment error, pro-viding users with more detailed statistical information that re-late to the quality of the fragment alignment. Moreover, thesystem makes it possible to visualize the location and magni-tude of the alignment errors as a spatial distribution on the frag-ment surface. Interactions are available from the histogram plotwhich allow the user to select a range of error values. Onceselected, the errors associated with these values are visualizedspatially across the surface of the fragment as shown in fig-ure 9. These tools are valuable for understanding the qualityof the geometric reconstruction results and understanding thecomplex geometric inter-dependencies of a highly-fragmentedbone fracture.
Figure 10: A severity assessment report is shown which includes values for sev-eral quantities that are known to be key factors in determining fracture severity.(a) (b)Figure 9: (a) shows a histogram of alignment errors where the bin of verticeshaving an error of approximately 0.25 mm has been selected (green). (b) showsa visualization of the spatial distribution surface points (blue spheres) whichhave these errors on the fragment surface. Severity assessment for bone fractures is heavily influencedby a number of related key factors. One contribution of thispaper is to provide quantitative values for some of these keyfactors which are heretofore unavailable in any other fractureanalysis software. Figure 10 shows an example of the severityassessment report which includes computed values of key fac-tors for each fragment in case one. The following list describesthe values provided in the severity assessment report in detail:1.
Fragment Volume:
This is 3D volume of the region en-closed by the fragment surface (shown as “Size” in theseverity report table). This information is useful in deter-mining if the fragment is structurally stable and useful inclinical reconstruction, i.e., can it be used for fixation, oris too small to be used in reconstruction.2.
Fragment Fracture Surface Area : The area of the frag-ment fracture surface is a major factor in determining frac-ture severity as shown in (Beardsley et al., 2005). The areaof the fracture surface generated is directly related to theenergy that the fractured bone absorbed during the fractureevent.3.
Fragment Displacement:
The fragment displacement isthe 3D translation vector that moves the centroid of thefragment, which indicates how far a fragment has movedor dispersed during the fracture event.4.
Angular Dislocation:
The angular dislocation of the frag-ment is the angle between the principal axis of the bonefragment in its original position and the principal axis ofthe bone fragment in its reconstructed position.These quantities are closely linked to fracture severity and thedisplayed values may allow users to more accurately and objec-tively estimate fracture severity. Computation of fracture sur-face areas from 2D or 3D CT images are di ffi cult and the results Computerized Medical Imaging and Graphics (2021) 11 are often unreliable. Here, the total fracture surface area canbe calculated more accurately after the 3D fragment surfacesare segmented and anatomically classified. Physicians have noway to quantitatively estimate fragment displacement and angu-lar dislocation from image data. Current approaches rely uponthe physician’s visual assessment and experience. For high en-ergy fracture cases, accurately assessment of these values aredi ffi cult. Since the proposed system has unique capabilities tocompute the original and reconstructed position for each bonefragment, the key factors shown are inaccessible from any othersource. Values uniquely available from the proposed softwareare fragment displacement, angular dislocation, and the surfacearea for the fragment outer surface and fracture surfaces.
4. Results
The bone reconstruction system was used to reconstruct sixclinical fracture cases which range from low energy fractureevents such as 1.5 foot fall, to high energy fracture events suchas a 50 mph car accident. The patient data, injury cause, and theOrthopedic Trauma Association (OTA) classification for eachcase is shown in Table 1. Since these are real clinical cases, thepatient names have been removed to protect their privacy. Allof the six cases were assigned a numerical severity score rang-ing from (1-100) by three well-trained surgeons based on theirpersonal experience and subjective inference shown in column C , C , and C of Table 1. Of particular note is the high vari-ance across clinicians of these classifications, consistent withthe findings of (Humphrey et al., 2005), highlighting the needfor objective measures of fracture severity.The section 4.1 shows the additional reconstruction resultsfor each case. For all six cases, the fragment and intact surfaceswere segmented, partitioned and classified using the methodsdescribed in 3.3, 3.4 and 3.5 to provide the data necessary topuzzle-solve each fracture case.Figure 11 shows the displaced positions and three di ff er-ent views of reconstructed fragments for all six clinical frac-tures. By visual assessment of the reconstruction results, allcases were reconstructed successfully with the exception ofcase three.Column 3 of Table 2 summarizes the global alignment er-rors for all six cases. From the table, we can see that globalalignment error after the construction for all six cases are rela-tively small ( < Figure 11: Six clinical tibia plafond fractures are puzzle-solved. The originalfractured positions for the fragments are shown in the left column and threedi ff erent views of reconstructed fragments are shown in the remaining columns.2 / Computerized Medical Imaging and Graphics (2021)
Case / OTA classification Injury mechanism C1 C2 C3 Avg1 F 38 C32 MVA (50 mph) 60 55 60 582 M 21 B13 Fall (30 ft) 50 60 58 563 F 42 C21 MVA (30 mph) 62 80 79 744 M 20 C13 ATV 6 15 32 185 M 24 C23 Fall (12 ft) 55 57 62 596 M 34 C11 Fall (18 ft) 70 65 77 71
Table 1: Patient data, injury cause, OTA classification, and severity scores by three surgeons for each case. The higher severity scores indicate higher fractureseverity.
Case ID Completiontime (sec) Number ofpoints onintact GlobalAlignmentError (mm)1 140 24935 0.232 220 45529 0.273 272 50539 0.324 90 33630 0.345 430 68160 0.336 650 117549 0.27
Table 2: Puzzle solving performance: the time needed to run puzzle-solvingalgorithm for the fracture case on a 2.4GHz, 4GB RAM laptop and the globalalignment error. recorded reconstruction time is influenced by several variablessuch as the number of spin images computed from the intacttemplate and the fragment surfaces, number of fragments, num-ber of iterations of the ICP algorithms, etc. In §4.2 we discussa small modification made to boost performance. Table 4 in theAppendix records the results of each fragment processed by thesystem.
Each of the six cases processed by our system supplies valu-able insight to the versatility of the reconstruction method. Allcases except Case 3 (which will be explained below) weredeemed successful due to good visual geometric agreement andvia the value of the global alignment error. Case 1 containedten unique fragments to place, yet still achieved a global er-ror of 0.23mm. This demonstrates the e ffi cacy of the proposedmethod even in complex fracture cases. Case 2’s fragment A4underwent a displacement magnitude of almost five centime-ters (48.98mm), which is substantially large compared to thetotal size of the tibia’s epiphysis. Despite this, Case 2’s re-construction was successful, demonstrating the resilience of theproposed reconstruction approach to local minima. Case 3 wasthe only unsuccessful reconstruction of all test cases. Its frag-ment A2 clipped through the base fragment (A1) despite theanatomical impossibility. This may be attributed to the low res-olution of Case 3’s CT DICOMs. It had the lowest resolutionat 0.683mm per pixel which would increase the likelihood offalse positive matches in the generated spin images. Case 4 and5 consist of a small number of fragments (3 and 4 respectively).This resulted in higher points per fragment to match, increasingthe computational cost of the puzzle solution. However, recon- struction for this case was successful and demonstrates the sys-tem’s capacity for processing large fragment surface matches.Case 6 was unique in its low percentages of outer surface areaon each fragment. The reconstruction was still successful de-spite the smaller amount of points to match per fragment. As mentioned in §3.6.11, the mean curvature histogram bi-ased search algorithm significantly improves the speed of theautomatic puzzle-solving algorithm by reducing the number ofspin images computed for both the intact template surface andthe fragment surfaces. The puzzle-solving algorithm is a com-plex process which consists of many steps, and the time spentfor the reconstruction is a ff ected by intermediate steps such ascomputing the spin images for the intact template and the frag-ment surfaces, matching spin images, and aligning fragmentsurfaces to the intact template. In order to better understandthe improvements, the following equation (4) details the timespent for reconstruction for each step. t tl = t o + n i t s + (cid:16) t s n f + t m + t f + t a (cid:17) M f (4)In this equation, t o denotes the time spent for processing theintact template before computing the spin images on the intacttemplate such as surface sampling, computing occupied regionsand computing the mean curvatures for the intact template. Theterm n i t s denotes the time spent to identifying feature points onthe intact template and compute their spin images. The term t s n f denotes the time spent computing descriptors on the frag-ment surfaces. The term t m denotes time spent generating hy-pothesized surface correspondences which is a ff ected by n f and n i . The term t f denotes the time spent removing false matches.The term t a denotes time spent aligning the fragment surface tothe intact template which is a ff ected by number of hypothesizedcorrespondences being tested. The term K denotes the numberof fragments in the fracture case. The major factors that im-pact the total reconstruction time are n f , the number of spinimages computed on the fragment surface, and n i , the numberspin images computed on the intact template. The mean cur-vature histogram approach in §3.6.11 reduces the total recon-struction time by reducing both n i and n f significantly Table3 shows the quantitative values for total reconstruction time, t tl , average time spent for matching spin images and filteringmatches, t m + t f , and n i for each case. Note that accelerationprovides roughly an order of magnitude decrease on compu-tation time for matching and reduces the total reconstructiontime, including user interaction, by about 50%. Experimental Computerized Medical Imaging and Graphics (2021) 13
Case t tl w / o accel.(min) t tt w / accel.(min) Avg matchtime w / oaccel.(sec) Avg matchingw / accel.(sec) n i w / o accel. n i w / accel. Table 3: Computational performance observed for each case using the system enhancements discussed in §3.6.11. reconstructions were conducted using a laptop computer with2.4GHz dual core CPU with 4GB memory.Our accelerated reconstruction times are significantly shorterthan those reported in similar systems. Work in (Fürnstahl et al.,2012) reported an average total reconstruction time of 83 min-utes which is approximately 5 times longer than our worst-caseclinical reconstruction (16 mins) and more than 10 times longerthan the average time it requires for reconstruction across our 6cases (~5 mins 45 secs). Unfortunately the performance in timeof the methods described in (Okada et al., 2009; Moghari andAbolmaesumi, 2008) are not reported.
5. Conclusion
The proposed system is capable of virtually reconstructingbroken bone fragments for complex bone fracture cases, whichis currently an unsolved problem in automatic puzzle-solvingalgorithms and di ffi cult to achieve using manual methods. Thebone reconstruction system designed in this paper enables usersto understand fracture cases from both 2D (CT image) and 3D(fragment surface) imagery. The system represents a uniquecombination of state-of-the-art 2D /
3D image processing andsurface processing algorithms. The software is a comprehen-sive reconstruction tool that guides users from the first step, i.e.,segmenting raw CT image data, to the last step, i.e., generatingquantitative, critical information about the fracture’s severity.Finally, 3D visualization of fragment surfaces can provide im-portant information for surgical treatment, especially for artic-ular fractures which often have a poor prognosis.The system demonstrates the e ffi cacy of spin image recon-struction in bone fractures from various fracture events with awide variety of traits. Reconstruction was successful in bonefractures with both many and few fragments (Cases 1, 4, and5), fractures where fragments have undergone large displace-ments (Case 2), and fractures where little intact outer surfacearea remains (Case 6). The reconstruction only su ff ers in caseswhere the source CT images were generated at low resolutions(greater than 0.5mm per pixel, Case 3). In each case, globalalignment errors less than 0.35mm were achieved.With the detailed fracture analysis o ff ered by this tool, itcould potentially improve surgical treatment, as previously ex-plored in (Thomas et al., 2010, 2011). The system represents asignificant step towards assisting physicians in classifying frac-ture severity which is especially di ffi cult in high-energy fracturecases. Use of quantitative metrics in conjunction with visual assessment promises to reduce variability in fracture severityclassification and may serve to build consensus in these di ffi -cult cases. The puzzle-solving algorithms and their integrationas a software system represent significant advancements towardimproving the treatment of comminuted tibial plafond fracturesand it is entirely possible this technology can be applied to otherproblematic limb fracture cases. While not directly discussedin this article we also note that the computational 3D puzzlesolving framework provides a heretofore unavailable patient-specific blueprint for fracture reconstruction planning. Havinga suitable blueprint for restoring the original anatomy, it be-comes possible for the surgeon to pre-operatively explore lessextensive surgical approaches and to attempt new minimally in-vasive surgical approaches. Appendix A. Summary of per-fragment results for recon-struction
Table 4 records the results of each fragment processed by thesystem. Displacement and displaced angle demonstrate howfar each fragment moved from its reconstructed position in thefracture event. The fracture surface area provides insight to theenergy of the event. The percentage of outer surface area to totalsurface area provides insight to the di ffi culty of finding matchesto the intact bone template. References
Albrecht, T., Vetter, T., 2012. Automatic fracture reduction,in: Lecture Notes in Computer Science. Springer BerlinHeidelberg, pp. 22–29.Anderson, D., Mosqueda, T., Thomas, T., 2008a. Quantify-ing tibial plafond fracture severity: Absorbed energy andfragment displacement agree with clinical rank ordering.Journal of Orthopaedic Research 26, 1046–1052.Anderson, D.D., Mosqueda, T., Thomas, T., Hermanson, E.L.,Brown, T.D., Marsh, J.L., 2008b. Quantifying tibial pla-fond fracture severity: Absorbed energy and fragment dis-placement agree with clinical rank ordering. Journal of Or-thopaedic Research 26, 1046–1052. doi: .Arthur, D., Vassilvitskii, S., 2006. Worst-case and smoothedanalysis of the ICP algorithm, with an application to thek-means method, in: 47th Annual IEEE Symposium onFoundations of Computer Science, IEEE. doi: .4 / Computerized Medical Imaging and Graphics (2021) F r a g m e n t D i s p l ace m e n t ( mm ) D i s p l ace d A ng l e ( d e g ) O u t e r S u rf ace A r ea ( mm ) F r ac t u r e S u rf ace A r ea ( mm ) % O u t e r t o T o t a l A r ea Case1 A . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . Case2 A . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . Case3 A . A . . . . . A . . . . . A . . . A . . . . . A . . . . . Case4 A . A . . . . . A . . . . . Case5 A . . . A . . . . . A . . . . . A . . . . . Case6 A . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . A . . . . . T a b l e : S u mm a r yo fr e s u lt s f o r a ll fr a g m e n t s o f t h e s i x t e s t ca s e s . Computerized Medical Imaging and Graphics (2021) 15Beardsley, C.L., Anderson, D.D., Marsh, J.L., Brown, T.D.,2005. Interfragmentary surface area as an index of com-minution severity in cortical bone impact. Journal ofOrthopaedic Research 23, 686–690. doi: .Besl, P.J. & McKay, N., 1992. A method for registration of3-d shapes. IEEE Trans. PAMI 14, 239–256.Brumback, R., Jones, A., 1994. Interobserver agreement in theclassification of open fractures of the tibia. The results of asurvey of two hundred and forty five orthopaedic surgeons.The Journal of Bone and Joint Surgery 76, 1162–1166.Charalambous, C., Tryfonidis, M., Alvi, F., Moran, M., Fang,C., Samaraji, R., Hirst, P., 2007. Inter- and intra-observervariation of the schatzker and AO OTA classifications oftibial plateau fractures and a proposal of a new classifica-tion system. The Journal of Bone and Joint Surgery 89,400–404.Chen, Y., Medioni, G., 1991. Object modeling by registrationof multiple range images, in: Proceedings. 1991 IEEE In-ternational Conference on Robotics and Automation, IEEEComput. Soc. Press. doi: .Chowdhury, A., Bhandarkar, S., Robinson, R., Yu, J., 2009.Virtual multi-fracture craniofacial reconstruction usingcomputer vision and graph matching. Computerized Med-ical Imaging and Graphics 33, 333–342.Cimerman, M., Kristan, A., 2007. Preoperative planning inpelvic and acetabular surgery: The value of advanced com-puterised planning modules. Injury 38, 442–449.van Eijnatten, M., van Dijk, R., Dobbe, J., Streekstra, G.,Koivisto, J., Wol ff , J., 2018. CT image segmentation meth-ods for bone used in medical additive manufacturing. Med-ical Engineering & Physics 51, 6–16. doi: .Felkel, P., Bruckschwaiger, M., Wegenkittl, R., 2001. Imple-mentation and complexity of the watershed-from-markersalgorithm computed as a minimal cost forest. ComputerGraphics Forum 20, 26–35. doi: .Felson, D.T., Zhang, Y., Hannan, M., Anderson, J., 1993. Ef-fects of weight and body mass index on bone mineral den-sity in men and women, the framingham study. Journal ofBone and Mineral Research 8, 567–573.Fisher, R.A., 1915. Frequency distribution of the values ofthe correlation coe ffi cient in samples from an indefinitelylarge population. Biometrika 10, 507–521. URL: .F.Meyer, 1994. Topographic distance and watershed lines.Signal Processing 38, 113–125.Freeman H, G.L., 1964. Apictorial jigsaw puzzles: The com-puter solution of a problem in pattern recognition. IEEETransactions on Electronic Computers EC-13, 118–127.Fürnstahl, P., Székely, G., Gerber, C., Hodler, J., Snedeker,J.G., Harders, M., 2012. Computer assisted reconstruc-tion of complex proximal humerus fractures for preop-erative planning. Medical Image Analysis 16, 704–720.doi: .Hansen, C., Johnson, C., 2005. The Visualization Handbook.Referex Engineering, Butterworth-Heinemann.Harders, M., Barlit, A., Gerber, C., 2007. An optimized surgi-cal planning environment for complex proximal humerusfractures. MICCAI Workshop on Interaction in MedicalImage Analysis and Visualiszation , 104–110.Humphrey, C.A., Dirschl, D.R., Ellis, T.J., 2005. Interob-server reliability of a ct-based fracture classification sys-tem. Journal of orthopaedic trauma 19, 616–622.Johnson, A., Hebert, M., 1996. Recognizing Objects byMatching Oriented Points. Technical Report CMU-RI-TR-96-04. Robotics Institute. Pittsburgh, PA.Johnson, A.E., Hebert, M., 1997. Surface registration bymatching oriented points, in: Proceedings of the Inter-national Conference on Recent Advances in 3-D DigitalImaging and Modeling, IEEE Computer Society, Washing- ton, DC, USA. pp. 121–128.Kovler, I., Joskowicz, L., Weil, Y.A., Khoury, A., Kron-man, A., Moshei ff , R., Liebergall, M., Salavarrieta, J.,2015. Haptic computer-assisted patient-specific preopera-tive planning for orthopedic fractures surgery. InternationalJournal of Computer Assisted Radiology and Surgery 10,1535–1546. doi: .Kronman, A., Joskowicz, L., 2013. Automatic bone fracturereduction by fracture contact surface identification and reg-istration, in: 2013 IEEE 10th International Symposium onBiomedical Imaging, IEEE. doi: .Liu, P., 2012. A system for computational analysis and recon-struction of 3D comminuted bone fractures. Ph.D. thesis.UNC Charlotte.Looker, A.C., Melton, L.J., Harris, T., Borrud, L., Shepherd,J., McGowan, J., 2009. Age, gender, and race / ethnic dif-ferences in total body and subregional bone density. Os-teoporosis International 20, 1141–1149.Lorensen, W.E., Cline, H.E., 1987. Marching cubes: a highresolution 3D surface construction algorithm. ComputerGraphics 21, 163–169.Low, K.L., 2004. Linear Least-Squares Optimization forPoint-to-Plane ICP Surface Registration. Technical ReportTR04-004. University of North Carolina at Chapel Hill.Raleigh, NC.Marsh, J., Buckwalter, J., Gelberman, R., 2002. Articular frac-tures: Does an anatomic reduction really change the result?The Journal of Bone Joint Surgery 84-A(7), 1259–71.Moghari, M.H., Abolmaesumi, P., 2008. Global registra-tion of multiple bone fragments using statistical atlas mod-els: Feasibility experiments, in: 2008 30th Annual Inter-national Conference of the IEEE Engineering in Medicineand Biology Society, IEEE. doi: .Okada, T., Iwasaki, Y., Koyama, T., Sugano, N., Chen, Y.W.,Yonenobu, K., Sato, Y., 2009. Computer-assisted preop-erative planning for reduction of proximal femoral frac-ture using 3-d-CT data. IEEE Transactions on Biomedi-cal Engineering 56, 749–759. doi: .Paulano, F., Jiménez, J.J., Pulido, R., 2014. 3d seg-mentation and labeling of fractured bone from CT im-ages. The Visual Computer 30, 939–948. doi: .Pulli, K., 1999. Multiview registration for large data sets,in: Second International Conference on 3-D Digital Imag-ing and Modeling (Cat. No.PR00062), IEEE Comput. Soc.doi: .Scheuering, M., Rezk-Salama, C., Eckstein, C., Hormann, K.,Greiner, G., 2001. Interactive repositioning of bone frac-ture segments, in: Proceedings of the Vision Modeling andVisualization Conference 2001, pp. 499–506.Shadid, W., Willis, A., 2013. Bone fragment segmentationfrom 3d CT imagery using the probabilistic watershedtransform, in: 2013 Proceedings of IEEE Southeastcon,IEEE. doi: .Shadid, W.G., Willis, A., 2018. Bone fragment segmentationfrom 3d ct imagery. Computerized Medical Imaging andGraphics 66, 14 – 27. doi: https://doi.org/10.1016/j.compmedimag.2018.02.001 .Thomas, T.P., Anderson, D.D., Marsh, L.J., Brown, T.D.,2007. Displaced soft tissue volume as a metric of commin-uted fracture severity, in: Annual Meeting of the AmericanSociety of Biomechanics, p. 309.Thomas, T.P., Anderson, D.D., Mosqueda, T.V., Hofwegen,C.J.V., Hillis, S.L., Marsh, J.L., Brown, T.D., 2010. Ob-jective CT-based metrics of articular fracture severity toassess risk for posttraumatic osteoarthritis. Journal ofOrthopaedic Trauma 24, 764–769. doi: .Thomas, T.P., Anderson, D.D., Willis, A.R., Liu, P., Marsh,6 / Computerized Medical Imaging and Graphics (2021)J.L., Brown, T.D., 2011. ASB clinical biomechanics awardpaper 2010. Clinical Biomechanics 26, 109–115. doi: .Willis, A., Anderson, D., Thomas, T., Brown, T., Marsh, J.L.,2007. 3D reconstruction of highly fragmented bone frac-tures, in: Proceedings of SPIE Medical Imaging, pp. 1–10.Willis, A., Cooper, D.B., 2008. Computational reconstructionof ancient artifacts. IEEE Signal Processing Magazine 25,65–83.Willis, A., Zhou, B., 2010. Ridge walking for 3D surfacesegmentation, in: Fifth International Symposium on 3DData Processing Visualization and Transmission (3DPVT),Paris, France.Wolfson, H., Schonberg, E., Kalvin, A., Lamdan, Y., 1988.Solving jigsaw puzzles by computer. Annals of OperationsResearch 12, 51–64. doi: .Zhang, Z., 1994. Iterative point matching for registrationof free-form curves and surfaces. International Jour-nal of Computer Vision 13, 119–152. doi:10.1007/BF01427149