Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

An Automatic Genetic Algorithm Framework for the Optimization of Three-dimensional Surgical Plans of Forearm Corrective Osteotomies

An Automatic Genetic Algorithm Framework for the Optimization of Three-dimensional Surgical Plans... Three-dimensional (3D) computer-assisted corrective osteotomy has become the state-of-the-art for surgical treatment of complex bone deformities. Despite available technologies, the automatic generation of clinically acceptable, ready-to-use preoperative planning solutions is currently not possible for such pathologies. Multiple contradicting and mutually dependent objectives have to be considered, as well as clinical and technical constraints, which generally require iterative manual adjustments. This leads to unnecessary surgeon efforts and unbearable clinical costs, hindering also the quality of patient treatment due to the reduced number of solutions that can be investigated in a clinically acceptable timeframe. In this paper, we propose an optimization framework for the generation of ready-to-use preoperative planning solutions in a fully automatic fashion. An automatic diagnostic assessment using patient-specific 3D models is performed for 3D malunion quantification and definition of the optimization parameters’ range. Afterward, clinical objectives are translated into the optimization module, and controlled through tailored fitness functions based on a weighted and multi-staged optimization approach. The optimization is based on a genetic algorithm capable of solving multi-objective optimization problems with non-linear constraints. The framework outputs a complete preoperative planning solution including position and orientation of the osteotomy plane, transformation to achieve the bone reduction, and position and orientation of the fixation plate and screws. A qualitative validation was performed on 36 consecutive cases of radius osteotomy where solutions generated by the optimization algorithm (OA) were compared against the gold standard solutions generated by experienced surgeons (Gold Standard; GS). Solutions were blinded and presented to 6 readers (4 surgeons, 2 planning engineers), who voted OA solutions to be better in 55% of the time. The quantitative evaluation was based on different error measurements, showing average improvements with respect to the GS from 20% for the reduction alignment and up to 106% for the position of the fixation screws. Notably, our algorithm was able to generate feasible clinical solutions which were not possible to obtain with the current state-of-the-art method. Keywords 3D Surgical Planning · Automatic · Forearm · Corrective Osteotomy · Multi-objective Optimization Research Highlights • Automatic diagnosis strategy based on bony landmarks. • Automatic placement of the fixation plate. • Two-stage weighted multi-objective optimization based on a genetic algorithm • Novel bone protrusion evaluation considering bone contact and surfaces gaps. • Patient-specific screw optimization based on bone density information. • Capability of considering all types of common osteotomies: single-cut, opening wedge, closing wedge. Abbreviations 2D: Two-dimensional 3D: Three-Dimensional CA: Computer-assisted CAD: Computer-aided design CASPA: Planning software developed at the Computer Assisted Research and Development Group CT: Computed tomography DoF: Degrees of freedom GPU: Graphics Processing Unit GS: Gold standard OA: Optimization algorithm SSM: Statistical shape models ROM: Range of motion 1. Introduction Post-traumatic healing in non-anatomical positions (malunions) or congenital deformations of bones can cause limitations in the range of motion (ROM) of the patient, generate pain and, if not treated properly, result in severe degenerative pathologies such as osteoarthritis (Nagy et al., 2008). The current gold standard for surgical treatment of these pathologies is the restoration of the normal anatomy in a surgical procedure known as a corrective osteotomy. In this procedure, the pathological bone is cut into two or more fragments, the fragments are realigned (clinical term: reduced) to their physiological position and stabilized with an osteosynthesis implant (Schweizer et al., 2010). However, the correction of bone malunions is highly patient-specific, requiring performing a complex 6-degree-of-freedom (DoF) correction (rotation and translation) for each bone fragment, in order to restore the physiological anatomy of the patient. Therefore, corrective osteotomies are technically challenging to perform without careful diagnosis and detailed preoperative planning of the procedure. Moreover, the successful outcome of a corrective osteotomy depends also on precise intraoperative navigation of the bone reduction (Fürnstahl et al., 2016). Conventional two-dimensional (2D) preoperative planning approaches based on X-ray and Computed Tomography (CT) images fail to correctly assess the inherently three-dimensional (3D) nature of bone malunions (Schweizer et al., 2010). Another drawback of 2D preoperative planning is that it cannot be used for surgical navigation, and surgeons must rely mainly on outdated, rudimentary surgical techniques to achieve the desired bone correction. Computer-assisted 3D preoperative planning addresses those problems. Several works have established 3D preoperative planning as the state-of-the-art technique for corrective osteotomies (Athwal et al., 2003; Dobbe et al., 2011; Fürnstahl, 2010; Fürnstahl et al., 2016), due to its clear benefits in patient treatment. It allows the quantification of malunions and its corrections in all 6 DoF. 3D preoperative planning offers also the possibility of precise pre-calculation of the osteotomy fixation using 3D representations of the fixation plates and fixation screws (Dobbe et al., 2011; Miyake et al., 2012a; Schweizer et al., 2010; Schweizer et al., 2013). Lastly, it enables the translation of the preoperative planning intraoperatively by means of surgical navigation, based either on patient-specific instruments (PSI) (Fürnstahl et al., 2016; Miyake et al., 2011; Murase et al., 2008; Omori et al., 2014) or on optical navigation systems (Andress et al., 2018). Thus, the introduction of 3D preoperative planning techniques marked a paradigm shift in patient treatment, allowing performing complex surgical procedures that would not be possible using conventional techniques (Athwal et al., 2003; Dobbe et al., 2014; Fürnstahl et al., 2016; Kunz et al., 2013; Roner et al., 2017; Schweizer et al., 2013; Schweizer et al., 2016; Zdravkovic and Bilic, 1990). In our institution, the current state-of-the-art preoperative planning of long-bone osteotomies encompass the following steps: in a first step, patient-specific 3D triangular surface models (hereinafter: 3D models) of the bones are generated based on the segmentation of the CT data of the patient (Fürnstahl et al., 2008), which is part of our standard clinical procedure (Fürnstahl et al., 2008). After obtaining the patient-specific 3D bone models, a diagnosis of the malunion is performed by comparison of the pathological bone model to a reconstruction target (Figure 1A). Usually, a mirrored model of the contralateral healthy side is used as the reconstruction target. The deformity analysis allows the definition of the needed osteotomy cuts along the pathological bone model, which are simulated as shown in Figure 1B. The resulting bone fragments can then be realigned into their anatomical position by realigning the bone fragments to fit the reconstruction target (Fürnstahl, 2010; Fürnstahl et al., 2016; Murase et al., 2008). The final step of the preoperative planning procedure is the simulation of the fixation of the osteotomy by integrating 3D models of the fixation plate and screws (Figure 1C). In total, the planning process for corrective osteotomies requires the definition of 18 DoF, without including the orientation and position of the screws. This preoperative plan can be translated into the operation room by means of patient-specific navigation instruments, designed according to the preoperative plan and later 3D-printed, allowing the surgeons to perform intraoperatively step-by-step the previously simulated procedure (Fürnstahl et al., 2016; Omori et al., 2014; Rosseels et al., 2018). Despite clear advantages of 3D preoperative planning, the generation of preoperative planning solutions requires the manual calculation of the aforementioned steps by trial and error, even when using dedicated 3D planning software (Fürnstahl et al., 2016; Roner et al., 2017; Schweizer et al., 2010). The development of a clinically feasible solution also requires close collaboration between surgeons, providing the clinical knowledge, and engineers, who have the technical expertise. As the availability of surgeons needed for consultancy is often very limited, the generation of an optimal preoperative solution for extra-articular long-bone osteotomies can add up to 4 hours (Fürnstahl et al., 2016) and might involve several iterations of manual adjustments. This incurs unwanted clinical costs as a consequence of the human workload spent on manual processes. Moreover, only a reduced number of clinically feasible solutions can be investigated due to the constrained clinical timeframe. Figure 1: Overview of state-of-the-art preoperative planning of long-bone osteotomies. (A) 3D bone models obtained through segmentation of CT data. The pathological bone model is shown in white and the reconstruction target in green. Diagnosis of the malunion is done by comparison of the two bones. (B) An osteotomy plane (cyan) is created and the pathological bone model is cut, generating the proximal (white) and distal (magenta) bone fragments. (C) Bone reduction is simulated by aligning the generated fragments to the reconstruction target, and the osteotomy is fixated with a fixation plate (in blue) and its corresponding screws (in gray). One possibility to overcome those challenges is the implementation of a computer-based planning approach, able to systematically and automatically generate optimal preoperative planning solutions. The automatic generation can significantly decrease clinical costs and reduce unnecessary interaction times of the collaborators. It could also improve the quality of patient treatment by considering a larger range of parameters and solutions in order to generate better preoperative planning solutions than those obtained by a human planner. Moreover, with the current trend of the health industry towards digitalization of patient data and treatment solutions, automatic methods become an essential asset for optimal clinical treatments. However, the implementation of an automatic optimization approach is a challenging task. Early on, Zdravkovic and Bilic (1990) proposed a computer-assisted preoperative planning framework for orthopedic surgeries in order to handle the complex tasks involved in the preoperative planning process. The latter, includes multiple nonlinear, discontinuous and non-differentiable planning objectives, making their mathematical manipulation difficult. Some objectives are contradicting and tightly associated with each other, which can cause the worsening of one objective while trying to improve another one. In the field of corrective osteotomies of long bones, only a few works have tried to tackle the automatic optimization problem (Belei et al., 2007; Carrillo et al., 2017; Schkommodau et al., 2005). Although these approaches have been promising and are pioneers in the filed of automatic optimization for surgical planning of orthopedic surgeries, they still lack of solutions that can be readily applied in clinical practice without further modifications. In this paper, we present a multi-staged, multi-objective optimization approach based on the artificial intelligent methods for the generation of ready-to-use preoperative planning solutions. The system is capable of calculating solutions considering all common bone malunions, osteotomy types, and surgical approaches. We have introduce the following key contributions with respect to our previous work (Carrillo et al., 2017):  An automatic diagnosis of bone deformation based on bony landmarks, which allows automatic narrowing of the optimization search space. In our previous work, the diagnosis of the bone deformation was assumed as given.  Automatic placement of the fixation plate using information provided by statistical shape models (SSM), in contrast to the manual definition of the feasible fixation areas used in (Carrillo et al., 2017).  A novel bone protrusion evaluation that supports the generation of better solutions by considering bone contact and surfaces gaps. The use of bone protrusion represents a more realistic clinical metric than the minimization of the cut surface used in our previous work.  In contrast to the heuristic approach employed in our previous work, in this paper, we have developed a patient-specific screw purchase optimization based on bone density information extracted from CT data.  The capability of considering all types of common osteotomies (single-cut, opening wedge, closing wedge) along the entire Radius, in order to generate solutions that would be difficult to achieve for a human planner within a reasonable time. The algorithm presented in (Carrillo et al., 2017) was only able to handle distal radius osteotomies and was not capable of generating single-cut solutions.  We have also reduced the computational effort of the strategy presented in (Carrillo et al., 2017) by applying a different multi-stage approach.  Finally, we have provided a profound clinical and mathematical explanation of all the optimization objectives, which was previously missing. The developed optimization framework was validated clinically against state-of-the-art preoperative planning solutions on a consecutive series of patients with forearm malunions, previously treated at our institution. In the following, we will give a brief overview of existing preoperative planning approaches for orthopedic surgeries. In section 2, we describe in detailed the proposed optimization framework. In section 3, dataset, experimental set-up, and results are presented. In section 4 a discussion about results, impact, and limitations of this work is given. Finally, in section 5, we draw final conclusions and give an outlook about future work. 1.1. Related Work The first step towards the generation of a complete 3D preoperative planning is the patient diagnosis done through the analysis of the bone malunion. In 3D preoperative planning, patient-specific 3D bone models generated from the segmentation of multi-planar data (Lorensen and Cline, 1987) are superimposed over a healthy reconstruction template using semi-automatic registration methods (Kawakami et al., 2002; Schenk et al., 2016; Schweizer et al., 2010; Vlachopoulos et al., 2017). Afterward, the bone deformation is quantified by means of clinical and mathematical measurements (Fürnstahl et al., 2016; Gosse et al., 1997; Murase et al., 2008; Subburaj et al., 2010), providing the basis for the preoperative planning of the surgical correction. Here, 4 main objectives have to be considered: osteotomy cut plane, reduction of the bone fragments, position of the fixation plate and direction and position of fixation screws. Current state-of-the-art of 3D preoperative planning for corrective osteotomies uses dedicated planning software for the manual calculation of each of these objectives. We refer the interested readers to (Fürnstahl, 2010; Fürnstahl et al., 2016; Schweizer et al., 2016; Vlachopoulos et al., 2015) for more information about the 3D planning tool. The latter facilitates the process of generation of a 3D preoperative plan, however the basic primitives operations needed for creating a preoperative plan are cutting Boolean operations and the possibility for interactive transformation, which would be also available in any dedicated CAD software. The process of creating each of the steps involved in the preoperative plan is difficult to control even by skilled engineers, as any change done to one of the parameters (e.g., position of fixation screws, inclination of osteotomy plane) must be manually propagated across the different objectives (Athwal et al., 2003; Belei et al., 2007; Bilic et al., 1994; Carrillo et al., 2017; Fürnstahl, 2010). Current-state-of-the-art planning (Fürnstahl et al., 2016; Roner et al., 2017; Vlachopoulos et al., 2015) is also incapable of handling contradicting objectives, meaning for example that an improvement in the accuracy of the bone reduction can subsequently deteriorate the position of the fixation plate or generate solutions with non-feasible osteotomy cuts. Similarly, an improvement in the position of the osteotomy cut can cause a less fitting position for the placement of the fixation plate. Existing automatic methods have only solved a reduced problem set, failing to include all objectives needed for a ready-to-use clinical solution , i.e., osteotomy cut plane, reduction of the bone fragments, position of the fixation plate and direction and position of fixation screws. Such is the case of previously described methods (Eck et al., 1990; Menetrey and Paul, 2004), where an automatic calculation was performed only for the osteotomy plane. Eck et al. (1990) calculated the osteotomy plane in femoral head reduction planning using a nonlinear optimization algorithm, based on a least-square-approximation solver. Menetrey and Paul (2004) did a similar parametrization of the osteotomy plane and wedge size for osteotomies around the knee, optimizing only the reduction alignment. Schkommodau et al. (2005) developed a multi-objective optimization strategy for corrective osteotomies of lower extremities that considered the following optimization objectives: leg length, translation, antetorsion, and angulation. Their method solved a simplified osteotomy sub-problem with only 12 DoF, which did not include the position of the fixation plate or the fixation screws into the optimization process. The multi-objective problem was solved by a sequential quadratic programming algorithm and considered the influence of all these objectives in the preoperative planning. The approach was later extended by Belei et al. (2007) to account for different osteotomy types (closing wedge, opening wedge, single cut) and to consider also double osteotomy solutions. These approaches represent the first attempt to automate the planning of corrective osteotomy. However, the planning was based on simplified geometry rather than on patient anatomy and it relied on intraoperatively calibrated fluoroscopic datasets. Learning-based methods have proven to be effective in similar applications of medical image processing techniques (Criminisi and Shotton, 2013; Esfandiari et al., 2018; Tschannen et al., 2016). In the field of shoulder arthroplasty, Tschannen et al. (2016) presented an automatic algorithm for preoperative planning of the resection plane for arthroplasty of the proximal humeral head, based on random regression forests. The approach allowed controlling the orientation, position, and size of the prosthetic humeral head in relation to the humeral shaft, using a direct mapping between the CT image and the parameters of the resection plane. The estimation of the plane position using CT data could be of interest to speed-up the convergence of optimization algorithms. Another interesting application is found in the field of spine surgery, where the problem of pedicle screw placement and pose estimation has been extensively studied (Farshad et al., 2017; Scheufler et al., 2011). Also, Esfandiari et al. (2018) proposed an algorithm based on convolutional neural networks for the 6-DoF estimation of the screw position and direction using intraoperative fluoroscopy data and estimation of the bone density. Similar approaches, taking into account the bone density information of the patient, should be considered in long-bone osteotomies to increase the quality of the osteotomy fixation. To the best of our knowledge, only our previous work (Carrillo et al., 2017) deals with an 18-DoF optimization problem for the generation of preoperative planning solutions of corrective osteotomies. The optimization approach presented in (Carrillo et al., 2017) has been taken as the basis for the core optimization algorithm (section 2.3.3) of the herein presented framework. 2. Methods We have developed an optimization framework for the generation of ready-to-use preoperative planning solutions for corrective osteotomies. A complete overview of our approach is given in Figure 2. The framework receives the 3D bone models of the patient, the reconstruction target and the fixation plate as an input (section 2.1). Afterward, an automatic diagnosis of the malunion is performed by identification of the pathological area, encoding also feasible plate regions and associated clinical constraints (section 2.2). This information is used by the optimization module, where the generation and optimization of preoperative planning solutions take place (i.e., the simulation of the osteotomy cut, realignment of the fragment and virtual fixation of the osteotomy) using a multi-stage optimization strategy (section 2.3). A genetic algorithm approach is used for optimization in which the objectives are encoded in a real-valued chromosome. The optimization is controlled by means of fitness evaluation functions, which provide quality measurements for the attainment of the different objectives. The output of the framework is a complete preoperative planning solution, including position and orientation of the osteotomy plane, the 6-DoF transformation to achieve the reduction, and the position and orientation of the fixation plate and screws. Our approach has been designed to cover all long-bone malunions, however, we have focused the implementation of our method on the radius due to the high rate of radius osteotomy procedures performed yearly at our institution. Figure 2: Overview of the automatic optimization framework for the 3D planning of long bone corrective osteotomy surgeries. 2.1. Input data generation In our institution, patient-specific 3D models of bone anatomy are used in the standard treatment workflow for complex malunions. For a better understanding, we briefly describe the model generation process although it is not an integral part of our optimization framework. Patient-specific 3D bone models of the forearm were generated from CT data (slice thickness 1 mm; 120 kV; Philips Brilliance 64 CT, Philips Healthcare, The Netherlands) using thresholding and region-growing algorithms of commercial segmentation software (Mimics Medical, Version 19.0, Materialise 2016, Leuven, Belgium). 3D models are generated using the Marching Cube algorithm (Lorensen and Cline, 1987) and given in the form of triangular surface meshes (stereolithographic models; STL) as described by (Roscoe, 1988). In our institution, the mirrored model of the contralateral bone is always used as the anatomical (healthy) reconstruction target. As the last input parameter, a feasible osteosynthesis implant (fixation plate) must be chosen, used for the fixation of the bone fragments after osteotomy. All 3D models were transformed into a common anatomical reference frame. In the case of the radius, we have used the anatomical coordinate system described in (Vlachopoulos et al., 2015). The coordinate system is � � � � ⃑ 𝐶𝐶𝐶𝐶 and it is oriented as shown in Figure 3A. Different anatomical and clinically relevant landmarks denoted by regions were defined and used throughout the entire pipeline. 7 point sets of landmarks regions, of 5 mm radius size, were automatically generated on the distal parts of both, pathological radius {𝐿𝐿 𝑃𝑃 , … ,𝐿𝐿 𝑃𝑃 } and reconstruction 1 7 target {𝐿𝐿 𝑅𝑅 , … ,𝐿𝐿 𝑅𝑅 } as depicted on Figure 3B. In our implementation, each set 𝐿𝐿 𝑃𝑃 and 𝐿𝐿 𝑅𝑅 contains a total of 1 7 𝑙𝑙 𝑙𝑙 𝑁𝑁 = 50 points, and each point set 𝐿𝐿 𝑃𝑃 is in correspondence to the reciprocal point set 𝐿𝐿 𝑅𝑅 . 𝑙𝑙 𝑟𝑟 � � � � � ⃗ Figure 3: Anatomical Axis and Bony Landmarks. (A) Anatomical coordinate system; coordinate axis 𝑪𝑪𝑪𝑪 is depicted by 𝒙𝒙 � � � � � ⃗ the red arrow and represents the radioulnar direction, the coordinate axis 𝑪𝑪𝑪𝑪 is depicted by the green arrow and represents 𝒚𝒚 � � � � � ⃗ the axial direction (long bone axis), and coordinate axis 𝑪𝑪𝑪𝑪 is shown in blue and indicates the dorsal-palmar direction. (B) 𝒛𝒛 Anatomical landmarks of the distal radius: (i) dorsal distal edge of the sigmoid notch; (ii) palmar distal edge of the sigmoid notch; (iii) center of distal radius styloid; (iv) center of scaphoid facet; (v) center of lunate facet; (vi) Lister’s tubercle and (vii) center of radius’ watershed line. 2.2. Automatic Diagnosis of Bone Malunion A common clinical practice for the 3D diagnosis of a bone malunion is to align the healthy proximal and distal joints of the pathological bone to the corresponding regions of the reconstruction target (Schweizer et al., 2010b). This process allows obtaining a visual approximation of the deformed area of the pathological bone and it defines also the most proximal and most distal regions of the bone at which an osteotomy cut can be performed. The bone malunion can then be quantified in 3D by the relative 4𝑥𝑥 4 transformation matrix between the proximally and distally aligned bone models (Figure 4). The transformation matrix encodes the translation � � � � ⃑ and rotational deviation of the malunion in all 6 DoF with respect to the anatomical coordinate system 𝐶𝐶𝐶𝐶 (Schweizer et al., 2010; Vlachopoulos et al., 2017). Based on this well-established principle, we have developed an automatic method for the calculation of the deformity region, the deformity profile, and the quantification of the bone malunion. Figure 4: Method for automatic diagnosis of malunions. (A-B) The pathological bone P is proximally (A) and distally (B) 𝒑𝒑 aligned to the reconstruction target R. (C) Example calculation of windows function for 𝑷𝑷 ; function 𝑫𝑫𝒑𝒑 is calculated along � � � � � ⃗ ( ) the 𝑪𝑪𝑪𝑪 of the bone, from 𝑳𝑳 to 𝑳𝑳 , with window function 𝑾𝑾 (∙) of thickness 𝒕𝒕 = 𝟏𝟏 𝒎𝒎𝒎𝒎 ,∞,∞ . (D) Deformity 𝒚𝒚 𝒙𝒙 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 curves; Lower envelopes 𝒍𝒍𝒍𝒍𝒍𝒍 (bold black signal) of 𝑫𝑫𝒑𝒑 (𝑷𝑷 ) (blue signal), and 𝒍𝒍𝒍𝒍𝒍𝒍 (bold magenta signal) of 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 𝑫𝑫𝒑𝒑 (𝑷𝑷 ) (gray signal), are calculated, corresponding to both proximal and distal alignments. (E) Deformity region is shown in orange, between the most proximal and most distal planes obtained from the deformity threshold 𝒅𝒅 applied to 𝑻𝑻𝑻𝑻 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 𝒍𝒍𝒍𝒍𝒍𝒍 and 𝒍𝒍𝒍𝒍𝒍𝒍 . In the automatic diagnosis, the pathological bone model 𝑷𝑷 is coarsely aligned to the reconstruction target 𝑹𝑹 using a landmark-based registration method, where each 𝐿𝐿 𝑃𝑃 of the pathological bone is registered to its 𝑙𝑙 𝒎𝒎𝒎𝒎 𝒎𝒎𝒎𝒎𝒎𝒎 reciprocal 𝐿𝐿 𝑅𝑅 of the reconstruction target, using standard iterative closest point (ICP) (Besl and McKay, 1992). 𝑙𝑙 Afterwards, the minimum and maximum points 𝐿𝐿 and 𝐿𝐿 are respectively calculated, on the bone length 𝑚𝑚 𝑚𝑚 � � � � ⃗ axis 𝐶𝐶𝐶𝐶 with respect to 𝑃𝑃 . Let 𝑦𝑦 𝛿𝛿 (𝑥𝑥 ,𝑌𝑌 ) = arg min� 𝑥𝑥 − 𝑦𝑦 � (Eq. 1) 𝑗𝑗 𝑦𝑦 ∈ 𝑌𝑌 𝑗𝑗 be the nearest neighbor function yielding the point 𝑦𝑦 of a point set or line 𝑌𝑌 with the smallest Euclidean distance 𝑗𝑗 ( ) to a point 𝑥𝑥 . 𝐿𝐿 and 𝐿𝐿 , and the line segment L 𝜆𝜆 between them can then be defined as 𝑚𝑚 𝑚𝑚 � � � � ⃗ � � � � ⃗ � � � � ⃗ � � � � ⃗ 𝐿𝐿 = 𝛿𝛿 � ,𝐶𝐶𝐶𝐶𝑝𝑝 � | 𝑝𝑝 = min� ∙ 𝐶𝐶𝐶𝐶 𝑝𝑝 � , 𝐿𝐿 = 𝛿𝛿 � ,𝐶𝐶𝐶𝐶𝑝𝑝 � | 𝑝𝑝 = max� ∙ 𝐶𝐶𝐶𝐶 𝑝𝑝 � , 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑝𝑝 ∈ 𝑃𝑃 𝑝𝑝 ∈ 𝑃𝑃 𝑖𝑖 𝑖𝑖 and 𝐿𝐿 (𝜆𝜆 ) = 𝐿𝐿 + 𝜆𝜆 (𝐿𝐿 −𝐿𝐿 ) | 0 ≤ 𝜆𝜆 ≤ 1, 𝑚𝑚 𝑚𝑚 𝑚𝑚 where · denotes the scalar product. Following the recommendation of (Vlachopoulos et al., 2017) for the size of the proximal and distal alignment regions, L(𝜆𝜆 ) is used to determine the points of 𝑃𝑃 that lie in the most proximal and distal 30% with respect to � � � � ⃗ 𝐶𝐶𝐶𝐶 : 𝑦𝑦 � � � � ⃗ � � � � ⃗ � � � � ⃗ { 𝑝𝑝 | 𝑝𝑝 ∙ 𝐶𝐶𝐶𝐶 ≤ 𝐿𝐿 (0.3)∙ 𝐶𝐶𝐶𝐶 } and { 𝑝𝑝 | ||𝑝𝑝 ∙ 𝐶𝐶𝐶𝐶 − 𝐿𝐿 || ≤ ||𝐿𝐿 (0.3) − 𝐿𝐿 ||} 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 These points are then used to align the proximal and distal joints of 𝑷𝑷 to 𝑹𝑹 (Figure 4A). The alignment is done using ICP registration (Besl and McKay, 1992) and the proximally and distally aligned models are denoted by 𝒑𝒑 𝒅𝒅 𝑷𝑷 and 𝑷𝑷 , respectively. 𝑝𝑝 𝑑𝑑𝑝𝑝 Afterward, the deformity profiles 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) (Figure 4D, blue signal) and 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) (Figure 4D, grey signal) are obtained from the deformity function 𝐷𝐷 𝑝𝑝 (𝑃𝑃𝐶𝐶 ) = �𝑊𝑊 (𝑙𝑙 (𝜆𝜆 ),𝑝𝑝 ) − 𝛿𝛿 � (𝑙𝑙 (𝜆𝜆 𝑊𝑊 ),𝑝𝑝 ),𝑊𝑊 (𝑙𝑙 (𝜆𝜆 ),𝑅𝑅 )� , (Eq. 2) 𝑚𝑚 𝑚𝑚 � | | 𝑃𝑃𝐶𝐶 | | 𝑚𝑚 =1,.., 𝑃𝑃𝑃𝑃 𝜆𝜆 =0,0.1,..,1 | | ( ) where 𝑃𝑃𝐶𝐶 is a pointset, 𝑃𝑃𝐶𝐶 is the cardinality of the pointset and 𝑊𝑊 · denotes the window function { | 𝑝𝑝 ∈ 𝑃𝑃𝐶𝐶 𝑐𝑐 −𝑡𝑡 < 𝑝𝑝 < 𝑐𝑐 +𝑡𝑡 } 𝑚𝑚 𝑚𝑚 ( ) 𝑊𝑊 𝑐𝑐 ,𝑃𝑃𝐶𝐶 = � . (Eq. 3) ∅ ∶ 𝑜𝑜𝑡𝑡ℎ 𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 𝑐𝑐 is the center of the window and 𝑡𝑡 is an adjustable user-defined thickness. As shown in Figure 4C, we have � � � � ⃗ used a thickness 𝑡𝑡 = 0.01 𝑚𝑚𝑚𝑚 , and we have taken 𝑡𝑡 and 𝑡𝑡 to be large enough to enclose 𝑷𝑷 and 𝑹𝑹 along 𝐶𝐶𝐶𝐶 𝑦𝑦 𝑚𝑚 𝑧𝑧 𝑚𝑚 � � � � ⃗ and 𝐶𝐶𝐶𝐶 . 𝑧𝑧 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 In order to reduce noise in the signals, we have implemented filtering using the lower envelopes 𝑜𝑜𝑒𝑒𝑙𝑙 (Figure 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑝𝑝 𝑑𝑑 4D, bold black line) and 𝑜𝑜𝑒𝑒𝑙𝑙 (Figure 4D, bold magenta line) of 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) and (𝑃𝑃 ) , respectively. Based on 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑜𝑜𝑒𝑒𝑙𝑙 and 𝑜𝑜𝑒𝑒𝑙𝑙 , a deformity threshold 𝑑𝑑 = 10% defines the part of the bone in which the deformity profile 𝑇𝑇 ℎ deviates more than 10% from their minimum value (Figure 4E). This region defines the most proximal and most distal sites for feasible osteotomy cuts. Finally, the malunion is quantified by a 4𝑥𝑥 4 homogenous transformation matrix 𝑀𝑀 representing the 𝑙𝑙 difference between the proximally and distally aligned 𝑷𝑷 . 𝑀𝑀 contains the information about the rotational and 𝑙𝑙 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 � � � � ⃗ � � � � ⃗ translational error of the bone malunion with respect to 𝐶𝐶𝐶𝐶 . The rotational deviation along each direction of 𝐶𝐶𝐶𝐶 is obtained by applying the Euler-Rodrigues formula (Schneider and Eberly, 2002) to the rotational component of 𝑀𝑀 , which yields the Euler angles 𝜑𝜑 ,𝜑𝜑 and 𝜑𝜑 . The translational deviation along each coordinate axis is 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑧𝑧 obtained from the translational parts 𝑇𝑇 ,𝑇𝑇 and 𝑇𝑇 of 𝑀𝑀 (Schneider and Eberly, 2002). Thus, the 6-DoF 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑙𝑙 � ⃗ quantification of bone malunion is given by the vector 𝑉𝑉 = � 𝜑𝜑 ,𝜑𝜑 ,𝜑𝜑 ,𝑇𝑇 ,𝑇𝑇 ,𝑇𝑇 � . 𝑚𝑚𝑚𝑚 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑚𝑚 𝑦𝑦 𝑧𝑧 2.3. Multi-stage optimization Our optimization approach comprises three parts. First, each clinical goal is characterized mathematically by means of optimization parameters, boundaries, and constraints. These parameters are used for the automatic generation of solutions at each optimization step, and the boundaries and constraints are used to control the search space of the algorithm. Our mathematical translation of the clinical goals into optimization objectives is explained in section 2.3.1. Secondly, the optimization algorithm requires a way to evaluate the performance of the optimization objectives. In our algorithm, this evaluation is performed by tailored fitness functions, developed for the minimization or maximization of the optimization parameters. The appropriate definition of the fitness function is crucial for the generation of feasible solutions and convergence of the algorithm. We give a detailed description of the fitness functions in section 2.3.2. Lastly, an optimization strategy must be chosen for the generation of solutions. In our case, the optimization algorithm had to be able to handle multiple objectives with 18 DoF and nonlinear constraints. We have implemented a multi-stage optimization approach based on the multi-objective genetic algorithm NSGA-II (Deb et al., 2000), further explained in section 2.3.3. 2.3.1. Optimization Objectives A complete preoperative solution for long bone osteotomies entails the calculation of 4 clinical objectives, as shown in Figure 5: (A) The position and orientation of the osteotomy cut plane (6 DoF), (B) the translation and rotation of the generated bone fragments for the reduction to the reconstruction target (6 DoF), (C) the position and orientation of the fixation plate (6 DoF) and (D) the position and orientation of the plate’s screws into the reduced fragments. The definition of each of these objectives requires more than a simple mathematical translation into adjustable parameters. We have to carefully decide on the parameters to control such that the parameter space can be minimal, and the controllability of the objectives can be directly evidence from a parameter change. In this section, we describe the strategy on the mathematical parameters, constraints, and boundaries associated with each objective. 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 Figure 5: Optimization objectives. (A) Position and orientation of the osteotomy cut plane (in cyan) along the pathological bone (in white). (B) Reduction alignment of the generated bone fragments (in white) with respect to the reconstruction target (in green). (C) Position of the fixation plate (in gray) with respect to the proximal (in white) and transformed distal fragment (in magenta). (D) Position and orientation of the screws of the fixation plate (in blue and gray), used for the stabilization and fixation of the bone fragments. A. Osteotomy Plane (Figure 5A). The osteotomy plane determines the size and location of the osteotomy and the orientation of the bone cut. It is represented by the position 𝑙𝑙𝑝𝑝 (3 DoF) and normal 𝑛𝑛� ⃗ (3 DoF), with a 𝑝𝑝𝑙𝑙 � � � � ⃗ total of 6 DoF. Regarding 𝑙𝑙𝑝𝑝 , we optimize only the position 𝑙𝑙𝑝𝑝 relative to the long bone axis 𝐶𝐶𝐶𝐶 (Table 1). 𝑦𝑦 𝑦𝑦 � � � � ⃗ In extra-articular osteotomies, only the translation of the osteotomy plane along the bone length axis 𝐶𝐶𝐶𝐶 has 𝑦𝑦 � � � � ⃗ � � � � ⃗ an effect on the osteotomy, thus the translation of the osteotomy plane with respect to 𝐶𝐶𝐶𝐶 and 𝐶𝐶𝐶𝐶 is 𝑚𝑚 𝑧𝑧 redundant. Together with the reduction alignment, the osteotomy plane determines the type of osteotomy and the shape of the wedge between the generated fragments. In our approach, we have considered all common osteotomy types (Figure 6): opening-wedge, closing wedge and single-cut osteotomy. Figure 6A shows the situation of the bone previous to realignment and after performing the osteotomy cut (in cyan). In an opening-wedge osteotomy (Figure 6B) a bone cut is made and the generated bone fragments are reduced to their anatomical position, yielding a wedge-shaped opening. A too large gap may prevent healing or may result in implant failure due to instability. To avoid this, large gaps are filled with structural bone graft in the surgery to support healing (Fernandez, 1982; Fürnstahl et al., 2016; Murase et al., 2008). Oppositely, a closing-wedge osteotomy (Figure 6C) refers to an osteotomy where a bone wedge has to be removed to complete the bone reduction. This may occur due to an overlap between the two generated bone fragments after their reduction to the reconstruction target (Fürnstahl et al., 2016; Murase et al., 2008). Our algorithm calculates first the necessary reduction of the bone fragments and afterward the overlap is defined by simulating a second planar cut. Additionally, our approach has also the feature of calculating single-cut osteotomies (Figure 6D) (Dobbe et al., 2017; Sangeorzan et al., 1989). A single cut osteotomy is always preferred against closing wedge solutions due to better bone contact and faster bone healing (Merle d’Aubigné and Descamps, 1952; Roner et al., 2017; Sangeorzan et al., 1989). This type of osteotomy is technically more complex to plan manually, as the bone reduction is achieved by sliding and rotating the generated bone fragments along the osteotomy cut plane (Sangeorzan et al., 1989). We achieve the support of single-cut osteotomies in our framework by constraining the closing wedge osteotomy planes to be parallel � � � � ⃗ and overlying. In this case, the wedge is limited to a maximum overlap of 1 mm along the 𝐶𝐶𝐶𝐶 . 𝑦𝑦 In each optimization stage, we use the polygon clipping algorithm of (Vatti, 1992) to cut model P by plane 𝒑𝒑𝒑𝒑 𝒅𝒅 𝒅𝒅𝒎𝒎 𝒕𝒕 (𝑙𝑙𝑝𝑝 ,𝑛𝑛 ), yielding potential candidates of the proximal and distal fragments 𝑷𝑷 and 𝑷𝑷 , respectively. 𝑝𝑝𝑙𝑙 The type of osteotomy is controlled by the maximum and minimum allowed wedge size using nonlinear constraints. Figure 6: Osteotomy Types. (A) Situation after performing the osteotomy cut (shown in cyan) and before reduction to the reconstruction target. The proximal fragment is shown in white and the distal fragment in magenta. (B) Opening wedge: distal fragment (magenta) has been reduced to the reconstruction target, generating an opening between both fragments. (C) Closing wedge: in order to complete the osteotomy, after simulated reduction of the distal fragment (magenta), a second planar cut must be done to remove the wedge of the overlapping bone region (shown in red). (D) Single cut: reduction of the distal fragment is achieved by fragment rotation and translation along the osteotomy plane. 𝒅𝒅 𝒅𝒅𝒎𝒎 𝒕𝒕 B. Reduction Alignment (Figure 5B). Once the osteotomy plane cut is executed, 𝑷𝑷 must be rotated and translated to match the given reconstruction target. This process of bone reduction represents one of the most important objectives for clinicians, because the precision of the reduction significantly influences the surgical outcome. Mathematically, the reduction alignment can be represented by a 4x4 homogenous transformation � ⃗ matrix (A ), constructed from the rotation 𝑅𝑅 = � ,𝜑𝜑 𝜑𝜑 ,𝜑𝜑 � and translation 𝑇𝑇 = � ,𝑇𝑇 𝑇𝑇 ,𝑇𝑇 � parameters f 𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑓𝑓 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 along each axis (6 DoF; Table 1). C. Position of Fixation Plate (Figure 5C). Another factor influencing the success of the treatment is the position of the fixation plate. A poorly fitting plate results in a less controlled procedure, decreases stability of the reduction and even leads to inaccurate reduction, delayed bone healing or failure of the fixation plate (Blecha et al., 2005; Fürnstahl et al., 2016). Mathematically, the position of the fixation plate along the bone entails a 6 DoF task, which can be encoded in the optimization by a 4x4 homogenous transformation matrix 𝒍𝒍𝒙𝒙 � ⃗ 𝐴𝐴 . The latter is constructed from the rotation 𝑅𝑅 = � ,𝜃𝜃 ,𝜃𝜃 𝜃𝜃 � and t he translation 𝑇𝑇 = � ,𝑇𝑇 𝑇𝑇 ,𝑇𝑇 � of 𝑝𝑝 𝑝𝑝 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑝𝑝 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 � � � � ⃗ the fixation plate (Table 1), relative to the reference coordinate system 𝐶𝐶𝐶𝐶 . D. Screw Position and Orientation (Figure 5D). The quality of the solution is also considerably dependent on the position and orientation of the screws, which are required to achieve a stable bone-plate interface. Several studies have shown how the length of the fixation screw can affect the stability of the osteotomy and how poorly placed screws could lead to soft-tissue inflammation and ligament tearing (Blecha et al., 2005; Dougherty et al., 2008; Reitman et al., 2004). Mathematically, the position and orientation of a fixation screw S is directly related to the position and orientation of the fixation plate, and can be described by � ⃗ S = h A � + R t � (Eq. 4) i s p s s i i i � ⃗ where h is the screw axis, A is the 4𝑥𝑥 4 homogenous transformation matrix of the fixation plate position, s p R is a 4𝑥𝑥 4 rotation matrix encoding orientation of the screw relative to the plate and t is a constant s s i i translation vector containing the relative position of the screw hole with respect to the plate (Table 1). For uniaxial locking screw plate systems, where the orientation of the screw is already given by the manufacturer, � � � � � ⃗ � � � � � ⃗ R is assumed to be constant, i.e., R = dır I , where dır is a constant direction vector. In the s s s 4x4 s i i i i � ⃗ optimization, the orientation of each screw is described by a vector θ , whose components are the Euler � ⃗ angles derived from R . The number of vectors θ is determined by the number of screws 𝑁𝑁 = |{𝑒𝑒 }| of the s s 𝑑𝑑𝑚𝑚 𝑚𝑚 i i fixation plate. Once the optimization objectives have been defined, the optimization parameters are encoded directly as decimal values into a parametrization vector. The vector is further denoted as chromosome 𝑥𝑥⃗, defined as � ⃗ ⃗ ⃗ 𝑥𝑥⃗ = � 𝑛𝑛 𝑙𝑙𝑝𝑝 𝑛𝑛 𝑛𝑛 𝜑𝜑 𝜑𝜑 𝜑𝜑 𝑇𝑇 𝑇𝑇 𝑇𝑇 𝜃𝜃 𝜃𝜃 𝜃𝜃 𝑇𝑇 𝑇𝑇 𝑇𝑇 θ 𝜃𝜃 … 𝜃𝜃 � , (Eq. 5) 𝑦𝑦 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 s 𝑑𝑑 𝑑𝑑 𝑧𝑧 𝑧𝑧 𝑥𝑥 𝑦𝑦 𝑧𝑧 1 2 𝑁𝑁 𝑠𝑠𝑖𝑖 which will be used by the fitness functions to control changes in the optimization objectives. Table 1: Summary of optimization objectives and their associated fitness functions, optimization parameters and constraints Optimization Fitness Optimization Challenges Parameters Constraints Objective Function Parameters Represented by 𝑨𝑨 ; 𝒇𝒇 Landmark-  Accuracy of joint  Rotation 𝑅𝑅 𝑓𝑓 based Reduction surface Maximal allowed Registration (𝜑𝜑 , 𝜑𝜑 , 𝜑𝜑 ) 6 𝑚𝑚 𝑦𝑦 𝑧𝑧 Alignment  Dependent on reduction error Error � � ⃗  Translation T landmark weighting (𝑓𝑓 ) (𝑇𝑇 , 𝑇𝑇 , 𝑇𝑇 ) 𝑓𝑓𝑚𝑚 𝑓𝑓𝑦𝑦 𝑓𝑓𝑧𝑧  Avoid longitudinal  For single-cut and intraarticular 𝑟𝑟𝑃𝑃 𝑃𝑃 𝑑𝑑𝑝𝑝𝑑𝑑𝑃𝑃  Position 𝑝𝑝𝑙𝑙  𝑛𝑛� ⃗ ∥ 𝑛𝑛� ⃗ 𝑝𝑝𝑙𝑙 𝑝𝑝𝑙𝑙 cuts (𝑝𝑝𝑙𝑙 ,𝑝𝑝𝑙𝑙 , 𝑝𝑝𝑙𝑙 ) Osteotomy Bone Protrusion 𝑚𝑚 𝑦𝑦 𝑧𝑧 6  Size of bone wedge  Osteotomy type and  Normal 𝑛𝑛� ⃗ Plane (𝑓𝑓 ) 𝑝𝑝 𝑙𝑙  Range given by location malunion (𝑛𝑛 , 𝑛𝑛 , 𝑛𝑛 ) diagnosis of section pl 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑥𝑥 𝑦𝑦 𝑧𝑧 influences reduction 2.2 alignment  Bone plate Represented by 𝑨𝑨 ; 𝒑𝒑  Gaps / steps between penetration Distance  Rotation 𝑅𝑅 𝑝𝑝 bone fragments and  Maximal allowed Fixation Plate - Position (𝜃𝜃 , 𝜃𝜃 , 𝜃𝜃 ) plate to be avoid 𝑚𝑚 𝑦𝑦 𝑧𝑧 6 distance to bone Fixation Plate Bone Fragments  Stable alignment � ⃗  Osteotomy site:  Translation 𝑇𝑇 𝑝𝑝 (𝑓𝑓 )  Various plate types palmar, dorsal, (𝑇𝑇 , 𝑇𝑇 , 𝑇𝑇 ) 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 lateral.  Proximity to joint � ⃗ ⃗ ⃗ � θ 𝜃𝜃 … 𝜃𝜃 � s 𝑑𝑑 𝑑𝑑  Minimal distance to 1 2 𝑁𝑁 𝑖𝑖𝑠𝑠 area 3 x 𝑁𝑁 𝑚𝑚𝑑𝑑 osteotomy plane Position and  Avoid penetration (𝑁𝑁 is the 𝑚𝑚𝑑𝑑 𝜃𝜃 = (𝜃𝜃 ,𝜃𝜃 ,𝜃𝜃 ) Screw Purchase 𝑑𝑑 𝑑𝑑 𝑑𝑑 𝑑𝑑  Number of screws 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑚𝑚 𝑦𝑦 𝑧𝑧 with osteotomy Orientation of number of (𝑓𝑓 ) inside the bone 4 for 𝑒𝑒 = 1, … ,𝑁𝑁 𝑚𝑚𝑑𝑑 Screws plane fixation  Bi-cortical solutions w.r.t plate coordinate  Patient-specific bone screws) preferred system and 𝐴𝐴 𝑝𝑝 density 2.3.2. Fitness functions and constraints Genetic algorithms employ minimization/maximization functions for the quantitative evaluation of the optimization parameters 𝑥𝑥⃗ during the optimization process. These functions are known as fitness functions and their definition is always challenging as they must be capable of providing a quantitative measurement of the quality of the different optimization objectives. The design of these fitness functions is a very challenging task, as we must consider not only mathematically solid functions that can be used by an optimization algorithm, but we must also guarantee that these functions reflect the clinical needs. For example, we cannot measure the quality of the osteotomy cut by simply measuring the inclination angle of the cut, because the bone contact and the bone protrusion of the generated fragments are of more clinical relevance than the orientation of the plane. In this section, we explain the design of the 4 fitness functions that we have developed based on clinically relevant measurements. Control of bone reduction alignment From the clinical point of view, some bone regions (i.e., joint regions) are more important to be precisely aligned to the reconstruction target than others. Moreover, in some cases, some bone regions inevitably deviate due to the pathology. We have developed a landmark-based error measurement to fine-control the reduction alignment by landmark regions 𝐿𝐿 𝑃𝑃 and 𝐿𝐿 𝑅𝑅 of the pathological bone and the reconstruction, respectively (Figure 3B). A 𝑙𝑙 𝑙𝑙 weight 𝑒𝑒 between 0 and 1 was assigned to all points of 𝐿𝐿 𝑃𝑃 , with ∑𝑒𝑒 = 1. Landmark regions located on the 𝑙𝑙 𝑙𝑙 𝑙𝑙 joint surface (i-v; Figure 3B) were given a bigger weight due to the importance of the joint for the reconstruction { } accuracy, which yields the following weight distribution 𝑒𝑒 , … ,𝑒𝑒 = 0.18 , 𝑒𝑒 ,𝑒𝑒 = 0.05 . In this way, the 1 5 6 7 quality of the reduction alignment is controlled by |{𝐿𝐿𝑃𝑃 }| 𝑁𝑁 =50 1 1 2 𝑓𝑓 = � 𝑒𝑒 � � � 𝑝𝑝 𝐴𝐴− 𝑞𝑞 � 𝑒𝑒 ℎ𝑒𝑒𝑒𝑒𝑒𝑒 𝑝𝑝 ∈ 𝐿𝐿 𝑃𝑃 ,𝑞𝑞 ∈ 𝐿𝐿 𝑅𝑅 (Eq. 6) 1 𝑙𝑙 𝑓𝑓 𝑚𝑚 𝑚𝑚 𝑚𝑚 𝑙𝑙 𝑚𝑚 𝑙𝑙 |{𝐿𝐿𝑃𝑃 }| 𝑁𝑁 𝑚𝑚 =1 𝑙𝑙 =1 which is the average weighted 𝑅𝑅𝐶𝐶𝑀𝑀𝑅𝑅 among the Euclidean distances between all 𝐴𝐴 -transformed points 𝑝𝑝 of the 𝑓𝑓 𝑚𝑚 pathological landmark regions and the points 𝑞𝑞 of the landmark regions of the reconstruction target. |{𝐿𝐿𝑃𝑃 }| is 𝑚𝑚 the cardinality of the set 𝐿𝐿𝑃𝑃 . The points are in correspondence, meaning that 𝑝𝑝 and 𝑞𝑞 have the same position 𝑚𝑚 𝑚𝑚 relative to their landmark centers. 𝐴𝐴 is constructed from 𝑅𝑅 = (𝜑𝜑 ,𝜑𝜑 ,𝜑𝜑 ) and 𝑡𝑡 = (𝑡𝑡 , 𝑡𝑡 ,𝑡𝑡 ) obtained 𝑓𝑓 𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑓𝑓 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 from chromosome 𝑥𝑥⃗. To guarantee a clinically acceptable reduction alignment, 𝑓𝑓 must be within a user-defined error boundary(≤ 0.5 𝑚𝑚𝑚𝑚 ). With this approach, it is possible to control the precision for specific regions on the bone. It allows not only a precise control of the reduction, but also the introduction of small deviations to the reconstruction target. Such deviations can support finding a better overall solution in case of contradicting targets by accepting small errors in the reduction alignment in favor of improving other objectives. For example, in case that the position of the fixation plate with respect to the 𝐴𝐴 -transformed bone fragment might not be 𝑓𝑓 within the acceptable solution range, the algorithm could allow a larger transformation error, e.g., 0.5 mm, and this small variation could already generate solutions with better fitting positions of the fixation plate. Control of osteotomy plane The parts of the bones that deviate from the anatomical target and that entail a surface gap among the fragments, are often referred to as bone protrusion (Figure 7A). The optimization of bone protrusion is very challenging because it is dependent on both, the osteotomy plane and the reduction alignment. We have decided to control the osteotomy cut position and direction implicitly by minimizing the bone protrusion among the generated fragments, in order to have a measurement that can be easily described to surgeons. As a positive side effect, the explicit optimization of the osteotomy plane is avoided, which has been proven to be challenging (Athwal et al., 2003; Schkommodau et al., 2005). The bone protrusion function 𝐵𝐵 𝑃𝑃 is defined as ( ) ( ) ( ) ( ) 𝐵𝐵 𝑃𝑃 𝐹𝐹 ,𝑅𝑅 = � 𝑊𝑊 𝑙𝑙𝑝𝑝 ,𝑝𝑝 − 𝛿𝛿 � 𝑙𝑙𝑝𝑝 ,𝑊𝑊𝑝𝑝 ,𝑊𝑊 𝑙𝑙𝑝𝑝 ,𝑅𝑅 � , (𝐄𝐄𝐄𝐄 .𝟕𝟕 ) 𝑚𝑚 𝑚𝑚 |𝐹𝐹 | 𝑝𝑝 ∈𝐹𝐹 𝑖𝑖 where 𝐹𝐹 is the point set of a bone fragment and 𝑅𝑅 is the reconstruction target. Windows thickness 𝑡𝑡 is set to � � � � ⃗ � � � � ⃗ 𝑡𝑡 = 4 𝑚𝑚𝑚𝑚 and 𝑡𝑡 and 𝑡𝑡 large enough to enclose 𝐹𝐹 and 𝑅𝑅 along 𝐶𝐶𝐶𝐶 and 𝐶𝐶𝐶𝐶 . 𝑦𝑦 𝑚𝑚 𝑧𝑧 𝑚𝑚 𝑧𝑧 In the optimization process, fitness function 𝑓𝑓 is defined as 𝑝𝑝 𝑟𝑟𝑝𝑝𝑚𝑚 𝑑𝑑 𝑚𝑚 𝐵𝐵 𝑃𝑃 (𝑃𝑃 ,𝑅𝑅 ) + 𝐵𝐵 𝑃𝑃 � 𝑃𝑃 𝐴𝐴 ,𝑅𝑅 � 𝑓𝑓 𝑓𝑓 = . (𝐄𝐄𝐄𝐄 .𝟖𝟖 ) 𝒑𝒑𝒑𝒑 An illustration of the bone protrusion calculation for 𝑷𝑷 is shown in Figure 7B-C. 𝒍𝒍𝒙𝒙 𝑑𝑑𝑑𝑑 Figure 7: Example of bone protrusion calculation. (A) Simplified drawing of the geometrical representation of bone protrusion after osteotomy cut and realignment. Gray cylinder represents the proximal fragment, orange cylinder represents 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 the reduced distal fragment. (B) Graphical explanation of the window function for fragment 𝑷𝑷 ; window function 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 ( ) 𝑾𝑾 𝒍𝒍𝒑𝒑 ,𝑷𝑷 is shown in blue, osteotomy plane is shown in cyan, the proximal fragment is shown in white and the 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 ( ) reconstruction target is shown in transparent green. (C) Calculation of the bone protrusion 𝑷𝑷𝑩𝑩 𝑷𝑷 ,𝑹𝑹 , which can be 𝐩𝐩 interpreted as the calculation of the RMSE between the point sets 𝐏𝐏 (white points) and 𝑹𝑹 (green points) subject to the function window 𝑾𝑾 (see Eq. 7). Control of plate position Stability of the bone reduction and, consequently, successful healing, depends on an appropriate implantation of the fixation plate. In general, a minimal distance between fixation plate FP and bone surface is desirable. However, in orthopedic surgeries, the position of the fixation plate is often constrained by the approach site of the surgery (Klausmeyer and Mudgal, 2014) and by the intrinsic anatomy of the bones (muscle and ligament attachments and joints). To integrate this information into our algorithm, we have defined clinically feasible plate location regions according to all standard approach sites (Klausmeyer and Mudgal, 2014). These regions were annotated by clinical experts on an average mean model, generated by a statistical shape model (SSM) of the forearm (Mauler et al., 2017; Sepasian et al., 2014) (Figure 8A). The annotation was performed using the Scalismo mesh painting tool (Graphics and Vision Research Group, University of Basel, Switzerland). Region transfer from the mean model to each patient-specific model was achieved through a model fitting registration 𝑚𝑚 algorithm described by (Lüthi et al., 2018). Let 𝐵𝐵 𝑅𝑅 ∈ 𝑃𝑃 the i-th feasible bone region fitted to the pathological ( ) bone that can be retrieved by a function 𝐵𝐵 𝑅𝑅 𝑃𝑃 ,𝑒𝑒 . 𝑚𝑚 As a first pre-processing step, the fixation plate is coarsely aligned to each 𝐵𝐵 𝑅𝑅 of 𝑃𝑃 through ICP methods (Besl and McKay, 1992) and used as input for the optimization algorithm. In a second preprocessing step, an that should be brought into contact automatic identification is performed to determine plate model points 𝐹𝐹𝑃𝑃 𝑃𝑃𝐹𝐹 � � � � ⃗ with the bone surface. To this end, an inherent coordinate system 𝐶𝐶𝐶𝐶 for the fixation plate is calculated using a principal component analysis (PCA; (Jolliffe, 2011)) of FP. The eigenvectors of the PCA correspond to the 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 � � � � ⃗ Euclidian coordinates of the 𝐶𝐶𝐶𝐶 and the origin 𝑐𝑐𝑒𝑒 is given by the mean of 𝐹𝐹𝑃𝑃 . Without loss of generality, let 𝑝𝑝 𝑃𝑃𝐹𝐹 � � � � ⃗ the positive direction of the 𝐶𝐶𝐶𝐶 axis point towards the undersurface of the plate. A point of 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 is 𝑧𝑧 𝑚𝑚 ∗ 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 � � � � ⃗ included in point set 𝐹𝐹𝑃𝑃 if (𝑝𝑝 − 𝑐𝑐𝑒𝑒 )∙ 𝐶𝐶𝐶𝐶 > 𝜀𝜀 𝑡𝑡 , where 𝑡𝑡 is the thickness of the fixation plate along 𝑚𝑚 𝑝𝑝 𝑧𝑧 𝑧𝑧 𝑧𝑧 𝑃𝑃𝐹𝐹 ∗ � � � � ⃗ 𝐶𝐶𝐶𝐶 . We empirically set 𝜀𝜀 = 0.85 , which worked for all plate types. An example of the identified 𝐹𝐹𝑃𝑃 point 𝑧𝑧 set is shown in Figure 8B. 𝐩𝐩𝐩𝐩𝐩𝐩 Figure 8: Feasible plate regions. (A) Feasible plate regions 𝑩𝑩 𝑹𝑹 (light yellow), encoded in the mean model of a radius 𝒎𝒎 SSM. The remaining colored regions are shown as indication of non-feasible plate regions (B) Surface points 𝑭𝑭𝑷𝑷 (in green) 𝑭𝑭𝑷𝑷 � � � � � ⃗ of the edge of the plate that should be in contact with the bone after registration. Coordinate axis 𝑪𝑪𝑪𝑪 is shown in black and 𝒛𝒛 is oriented towards the undersurface of the plate. 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 On each optimization step, the average distance between the fixation plate and the proximal fragment 𝑃𝑃 is calculated as 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 ( ) 𝐷𝐷 = �� 𝛿𝛿 � 𝐴𝐴 𝑝𝑝 ,𝑅𝑅𝐵𝐵 𝑃𝑃 ,𝑒𝑒 � − 𝑝𝑝 � , 𝑒𝑒 ℎ 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 (Eq. 9) ∗ 𝑝𝑝 𝑗𝑗 𝑗𝑗 𝑗𝑗 � 𝐹𝐹𝑃𝑃 � 𝑚𝑚 =1,..,|{𝐵𝐵 𝑅𝑅 }| 𝑖𝑖 𝑗𝑗 =1,…,� 𝑃𝑃𝐹𝐹 � ∗ ∗ where |𝐹𝐹𝑃𝑃 | is the cardinality of 𝐹𝐹𝑃𝑃 , 𝐴𝐴 is the 4x4 homogenous transformation matrix of the plate obtained 𝑝𝑝 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 from chromosome 𝑥𝑥⃗, 𝐵𝐵 𝑅𝑅 (𝑃𝑃 ,𝑒𝑒 ) are the points of feasible region 𝑒𝑒 contained in 𝑃𝑃 , and 𝛿𝛿 is the already defined nearest neighbor function of Eq. 3. Similarly, the distance between the fixation plate 𝐹𝐹𝑃𝑃 and the distal 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 fragment 𝑃𝑃 is obtained by 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝐷𝐷 = �� 𝛿𝛿 � 𝐴𝐴 𝑝𝑝 ,𝐴𝐴 𝑅𝑅𝐵𝐵 (𝑃𝑃 ,𝑒𝑒 )� − 𝑝𝑝 � , where 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 (Eq. 10) ∗ 𝑝𝑝 𝑗𝑗 𝑓𝑓 𝑗𝑗 𝑗𝑗 � 𝐹𝐹𝑃𝑃 � |{ }| 𝑚𝑚 =1,.., 𝐵𝐵 𝑅𝑅 𝑖𝑖 𝑗𝑗 =1,…,� 𝑃𝑃𝐹𝐹 � where 𝐴𝐴 is the 4x4 homogenous transformation matrix of the distal fragment obtained from chromosome 𝑥𝑥⃗. The 𝑓𝑓 average of both distance values is used for controlling the position and orientation of the plate, encoded into fitness function 𝑓𝑓 as 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑓𝑓 = � 𝐷𝐷 +𝐷𝐷 � ⁄2 (Eq. 11) Control of screw position and orientation The bone density of the patient bone considerably influences screw stability, as cortical and trabecular bone distribution differs along the bone (Kubiak et al., 2006; Reitman et al., 2004). The lengths of the screws play another important factor, as too short or too long screws might cause complications or injuries (Blecha et al., 2005; Spahn, 2004; Wall et al., 2012). Also, the position and direction of screws with respect to the bone density considerably affect the stability of the fixation (Dougherty et al., 2008; Weninger et al., 2010). Therefore, the bone density distribution will be considered for defining the screw positions. We have developed a method for generating a patient-specific bone density mask based on the original CT image 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 𝐼𝐼 containing the patient-specific normalized Hounsfield values (𝐼𝐼 ) (Figure 9A-B). Let Φ (𝑝𝑝 ) be a 𝐼𝐼 𝑚𝑚 transformation function transforming the 3D coordinates of a point 𝑝𝑝 to image space coordinates � ,𝐼𝐼 ,𝐼𝐼𝐼𝐼 � with 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑧𝑧 respect to an image 𝐼𝐼 . The outer cortical surface of the pathological bone model 𝑷𝑷 (orange overlay in Figure 9B) is 𝑃𝑃𝑇𝑇𝐿𝐿 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 rasterized into an empty image 𝐼𝐼 of same size as 𝐼𝐼 using the method described by (Pineda, 1988). A grey 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 value of 1000 𝑚𝑚𝑎𝑎𝑥𝑥 (𝐼𝐼 (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 )) is assigned to each rasterized voxel. Next, a dilation of size 1 pixel is (𝑚𝑚 ,𝑦𝑦 ,𝑧𝑧 ) 𝑃𝑃𝑇𝑇𝐿𝐿 applied to 𝐼𝐼 using a standard dilation kernel and the weighted bone-density mask shown in Figure 9C is defined 𝑒𝑒𝑒𝑒𝑒𝑒 𝐿𝐿 𝑊𝑊 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 ( ) ( ) ( ) 𝑀𝑀 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 = max� 𝑀𝑀 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ,𝐼𝐼 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 � , (Eq. 12) 𝑃𝑃𝑇𝑇𝐿𝐿 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 which returns the values of 𝐼𝐼 for (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ) lying on the cortical layer, values of 𝑀𝑀 (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ) lying inside the pathological bone (trabecular bone) and 0 otherwise. The masking process was performed in Matlab (MATLAB, 2017) using the Iso2Mesh toolbox (Fang and Boas, 2009). In the objective evaluation of the optimization process, penetration points 𝑒𝑒𝑛𝑛 and 𝑜𝑜𝑜𝑜 𝑡𝑡 (Figure 9D-E, cyan 𝑑𝑑 𝑑𝑑 𝑖𝑖 𝑖𝑖 points) are calculated between each screw 𝑒𝑒 and the bone fragments in reduced position. The penetration points 𝑚𝑚 � ⃗ are found by casting a ray formed by screw center A t (Eq. 4, with 𝑅𝑅 = I ) and screw axis h and p s 𝑑𝑑 4x4 s i 𝑖𝑖 i calculating the intersection with the bone as described in (Mount and Arya, 1998; Schneider and Eberly, 2002) . Afterwards, a sampling between 𝑒𝑒𝑛𝑛 and 𝑜𝑜𝑜𝑜 𝑡𝑡 is performed in 3D space along the screw axis using line equation 𝑑𝑑 𝑑𝑑 𝑖𝑖 𝑖𝑖 𝐿𝐿 � � =λ 𝑒𝑒𝑛𝑛 +λ � 𝑡𝑡 −𝑜𝑜𝑜𝑜 𝑒𝑒𝑛𝑛 � | λ = {0,0.1, … ,1} (Figure 9D, cyan and red spheres). The sampling 𝑃𝑃 𝑃𝑃 𝑑𝑑 𝑃𝑃 𝑑𝑑 𝑑𝑑 𝑃𝑃 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 approach in 3D space is computationally more efficient than evaluating the entire density in the image space. The average weighted density value for each screw 𝑒𝑒 is obtained by 𝑚𝑚 𝑊𝑊 𝐷𝐷 𝑉𝑉 = � 𝑀𝑀 � Φ 𝑊𝑊 � 𝑙𝑙 � λ � � . (E� q. 13) 𝑃𝑃 𝐶𝐶 𝐶𝐶 𝑀𝑀 𝑖𝑖 𝑒𝑒 𝑒𝑒 � λ � 𝐶𝐶 𝑒𝑒 λ ={0,0.1,…,1} 𝐶𝐶 𝑒𝑒 The position and direction of the screws is controlled by the average of all the weighted density values 𝑁𝑁 𝑆𝑆 𝑖𝑖 𝑓𝑓 = � 𝐷𝐷 𝑉𝑉 (Eq. 14) 4 𝑃𝑃 𝑖𝑖 𝑁𝑁 𝑃𝑃 𝑖𝑖 𝑚𝑚 =1 𝐶𝐶𝑇𝑇 Figure 9: Fitness calculation for screw position and direction using patient-specific bone density information. (A) 3D 𝑶𝑶𝒑𝒑𝒎𝒎𝑶𝑶 models of the pathological radius (in orange) and ulna (in white); grey plane indicates a CT slice. (B) Slice of 𝑰𝑰 in coronal view. The radius is indicated by the semi-transparent orange mask and the outer cortical layer by the orange line. (C) 𝑾𝑾 ( ) Resulting weighted bone-density mask 𝑴𝑴 𝒙𝒙 ,𝒚𝒚 ,𝒛𝒛 . (D) Sampling calculation of insertion points for Screw 𝑪𝑪 ; the cyan 𝒎𝒎 spheres represent 𝒎𝒎 𝒎𝒎 and 𝒍𝒍 𝒐𝒐𝒕𝒕 and the red spheres represent the sampled points along the screw axis. (E) Coronal view of 𝒅𝒅 𝒅𝒅 𝒎𝒎 𝒎𝒎 𝑾𝑾 𝑴𝑴 (𝒙𝒙 ,𝒚𝒚 ,𝒛𝒛 ) showing the corresponding points of 𝑳𝑳 (𝛌𝛌 ) in image space. 𝑪𝑪 𝒎𝒎 2.3.3. Core Algorithm The core optimization algorithm is a modified version of the multi-stage multi-objective genetic optimization approach previously described in our work (Carrillo et al., 2017), capable of handling multiple contradicting objectives and nonlinear constraint. Once each planning objective has been parameterized using the real-valued chromosome 𝑥𝑥⃗, an optimal Pareto set 𝑌𝑌 of planning solutions can be found by 𝑚𝑚 𝑚𝑚 𝑚𝑚 { ∶ 𝑦𝑦 = min ( 𝑓𝑓 ( ) ( ) ( ) ) } : ℝ → ℝ , (Eq. 15) 𝑌𝑌 = 𝑦𝑦 𝜖𝜖 ℝ 𝑥𝑥⃗ , 𝑓𝑓 𝑥𝑥⃗ , … , 𝑓𝑓 𝑥𝑥⃗ | 𝑓𝑓 𝑚𝑚 ∈ Χ 1 2 𝑚𝑚 𝑚𝑚 𝑚𝑚 𝑚𝑚 where ℝ is the solution space generated by the fitness function 𝑓𝑓 (𝑥𝑥⃗), ℝ being the parameter space and Χ the 𝑚𝑚 𝑚𝑚 |𝑚𝑚⃗| 18 set of optimal chromosomes generated by the algorithm. In our case, ℝ = ℝ ≥ ℝ depending on the 𝑚𝑚 4 number of screws for each plate, and ℝ = ℝ where 𝑚𝑚 is the number of fitness functions. However, standard NSGA-II (Deb et al., 2002) is only able to find symmetrically optimized solutions with 𝑚𝑚 equally weighted optimization objectives on the solution space ℝ . In our optimization problem, each objective must have a different importance, according to its clinical relevance. Based on (Carrillo et al., 2017), we have applied the following weighting schema: ⁄ 𝑖𝑖 −k 𝑚𝑚⃗ 𝑃𝑃 𝑖𝑖 𝑒𝑒 (𝑥𝑥⃗) = ∑ 𝑒𝑒 . (Eq. 16) 𝑚𝑚 The weighting function 𝑒𝑒 increases the sparsity of solutions 𝑌𝑌 around the utopia point (Miettinen, 1999). The sparsity of 𝑓𝑓 (𝑥𝑥⃗) along the Pareto set is controlled by the constant 𝑘𝑘 , which represent the complementary 𝑚𝑚 𝑚𝑚 percentage of the desired sparsity. For instance, if 𝑘𝑘 = 30, 100− 𝑘𝑘 = 70% of solutions are to be found closer 𝑚𝑚 𝑚𝑚 to the utopia point. We have defined, together with surgeons, the following optimal weighting schema: Reduction alignment as the objective with the highest priority, followed by the osteotomy plane, the plate position and the screw purchase, respectively. An initialization and automatic verification of constraints and boundaries is done previous to each optimization stage. The distribution of the optimization process into different stages allows handling contradictive objectives by giving different importance on each stage to different objectives. Initialization Stage 1. An initial population 𝑋𝑋 , where the superscript denotes the optimization stage, of 200 chromosomes is randomly generated among the parameter range of chromosome 𝑥𝑥⃗ (Eq. 5). At each stage, we optimize over all parameters (all objectives) but considering only a subset of fitness functions. A summary of the algorithm parameters is given in Table 2. Stage 1. In this stage, the reduction alignment (𝑓𝑓 ), the osteotomy plane (𝑓𝑓 ), and the position of fixation plate 1 2 (𝑓𝑓 ) are optimized. Our algorithm is run with 3 objectives and weighted towards 𝑓𝑓 (𝑘𝑘 = 20, Eq. 16) and 𝑓𝑓 3 1 1 2 (𝑘𝑘 = 40, Eq. 16). This strategy permits to obtain solutions for the reduction alignment and bone protrusion close to the utopia point, but allowing a larger freedom on the plate position (𝑘𝑘 is set to 100, i.e., sparsity is unaffected). By forcing the algorithm to remain in the solution space where the reduction alignment (𝑓𝑓 ) and the osteotomy plane (𝑓𝑓 ) are close to the utopia point, contradictive solution that might engender the quality of 𝑓𝑓 or 2 1 𝑓𝑓 , are automatically neglected by the evolutionary process of the GA. The resulting Pareto set is stored in matrix 1 1∗ 𝑌𝑌 containing the best 70 solutions corresponding to the Pareto front. Accordingly, 𝑋𝑋 represents the corresponding set of parameters yielding 𝑌𝑌 . Table 2: Summary of algorithm parameters for each stage Parameter Stage 1 Stage 2 Input Population 𝑋𝑋 𝑋𝑋 Max # Generations 200 200 Population Size 200 200 Number of objectives 3 3 𝑓𝑓 , 𝑓𝑓 , 𝑓𝑓 𝑓𝑓 , 𝑓𝑓 , 𝑓𝑓 Fitness Functions 𝑓𝑓 (𝑥𝑥⃗) 1 2 3 1 3 4 𝑚𝑚 𝐾𝐾 = 20, 𝑘𝑘 = 40 𝐾𝐾 = 40 Weighting Schema(𝐾𝐾 ) 1 2 3 𝑚𝑚 1 1∗ 2 2∗ Output 𝑌𝑌 , 𝑋𝑋 𝑌𝑌 , 𝑋𝑋 Initialization Stage 2. The best 70 solutions from stage 1 are used for the initialization of the first generation of stage 2. To this end, the parameter range 𝑋𝑋 of the initial population of stage 2 are randomly initialized per component, such that each component value is constrained to lie between the maximum and minimum values 1∗ 1∗ 1∗ obtained from 𝑋𝑋 . Thus, 𝑥𝑥⃗ ∈ � 𝑚𝑚𝑒𝑒𝑛𝑛 �� ,𝑋𝑋𝑚𝑚𝑎𝑎𝑥𝑥 �� 𝑋𝑋 wher� e 𝑒𝑒 is the component i-th of the chromosome 𝑥𝑥⃗. 𝑚𝑚 𝑚𝑚 𝑚𝑚 The algorithm parameters of stage 2 are given in Table 2. Stage 2. In this stage, we optimize the reduction alignment (𝑓𝑓 ), the position of fixation plate (𝑓𝑓 ) and the 1 3 average screw purchase (𝑓𝑓 ). Our algorithm is weighted towards 𝑓𝑓 (𝑘𝑘 = 40, 𝑬𝑬𝑬𝑬 .𝟏𝟏𝟏𝟏) , to guarantee solutions 4 3 3 with an optimal fixation plate alignment, and to avoid deteriorating the reduction alignment. The remaining fitness functions are not weighted (𝑘𝑘 = 𝑘𝑘 = 100). The resulting Pareto set is stored in matrix 𝑌𝑌 , containing 1 2 2 the best 70 optimal chromosomes corresponding to the Pareto front. Thanks to our weighing schema, 𝑌𝑌 contains 2∗ 2 optimal solutions that are all within an acceptable clinical errors. 𝑋𝑋 was ranked according to the solutions 𝑌𝑌 2∗ 𝑃𝑃𝑝𝑝𝑙𝑙 by the best combined fitness among the 3 objectives of stage 2, yielding 𝑋𝑋𝑅𝑅 . The final output 𝑥𝑥⃗ of the algorithm is obtained by taking the solution with the best plate alignment among the ranked list � 𝑓𝑓 (𝑥𝑥� � ⃗)− 𝑚𝑚𝑒𝑒𝑛𝑛 (𝑓𝑓 (𝑥𝑥� � ⃗))� 3 3 𝑃𝑃𝑝𝑝𝑙𝑙 𝑥𝑥⃗ = min . 𝑥𝑥� � ⃗ ∈ 𝑅𝑅𝑋𝑋 (� � ⃗) (� � ⃗) � 𝑚𝑚𝑎𝑎𝑥𝑥 (𝑓𝑓 𝑥𝑥 − 𝑚𝑚𝑒𝑒𝑛𝑛 (𝑓𝑓 𝑥𝑥 � 3 3 3. Results 3.1. Datasets We have performed a qualitative validation (section 3.2) and a quantitative evaluation (section 3.3) on retrospective cases of malunited radii, which have been included in a large clinical trial about CA corrective osteotomy. From these consecutive cases, 36 cases were eligible according to the inclusion criteria given in Table 3. All 36 patients were treated at our orthopedics department between 2015 and 2017 and underwent navigated forearm osteotomy surgery through 3D preoperative planning and patient-specific instrumentation. Informed consent to use the patient data in an anonymized form for computer analyses was obtained, with ethical approvals KEK-ZH-Nr. 2015-0186 and BASEC-Nr. 2018-01608. The corresponding CT scans used for the generation of the 3D models were obtained according to a standard scanning protocol (slice thickness 1 mm; 120 kV; Philips Brilliance 64 CT, Philips Healthcare, The Netherlands). 3D preoperative planning of all patients were prospectively created in a manual fashion using a preoperative planning software (CASPA, Balgrist CARD AG, Switzerland) by the responsible hand surgeon, together with an expert planning engineer. These solutions were considered as the Gold Standard (GS). The baseline of the selected data is presented in Table 3. Table 3: Inclusion and exclusion criteria for clinical and technical evaluation of the clinical cases Inclusion Exclusion  Signed informed consent for data use  Incomplete CT data  Revision surgery  Age ≥ 16 years old  Radius osteotomy  Intra-articular corrective osteotomy  Availability of both, pathological and  Multiple osteotomies on one bone contralateral CT data  Complete 3D preoperative planning including fixation plate Table 4: Baseline data of the 36 cases included in the study Mean SD Age (years)* 33 ±14 Male Female Gender 11 25 Right Left Affected Side 19 17 Distal Radius Radius Shaft Location Malunion 14 22 Single Closing Opening Cut Wedge Wedge Osteotomy Type 11 8 17 Fixation Plate 9 different types * years at time of surgery 3.2. Qualitative Validation A preoperative planning solution (optimization solution; OA) was generated by our algorithm for each of the 36 consecutive cases in an automated fashion. The anatomical site and osteosynthesis implant were defined to be the same as in the GS solution for each case. A clinical survey was performed with 6 readers (2 senior hand surgeons, 1 board-qualified orthopedic surgeon, 1 resident surgeon and 2 expert planning engineers). In the survey, we asked each reader the question: “which of the two presented solutions would you implement in the surgery without further modification?”. The two solutions (GS and OA) were presented to the readers in a blinded fashion and in random order. The assessment for each case is shown in Figure 10A. OA solutions were chosen as the better solution in 55% of the times (Figure 10B). Out of the 55%, surgeons’ votes represented 38% and engineers’ votes 17%. Similarly, for the 45% of the GS votes, 29% account for the votes of the surgeons and the rest 16% represent the votes of the engineers. Figure 10: Results of the qualitative validation. (A) Voting of each reader and case. Each case was evaluated by 4 surgeons (readers 1-4) and 2 experienced planning engineers (readers 5-6). (B) Percentage distribution among the different voters for GS and OA solutions. In Figure 11 we present two examples where the OA solutions outperformed the GS solutions. In both cases, the reduction alignment and the position of the fixation plate obtained by the OA were clearly better than the GS solutions. Moreover, in both cases, the distance of the fixation plate with respect to the bone fragments has been decreased and undesired bone intersections have been solved by the OA. In these two cases, 100% of the readers favored the OA solutions. There were also cases where the GS solutions were preferred by the majority of the readers. Figure 12 shows two example cases where more than 50% of the readers favored the GS solutions. In the case shown in Figure 12A, surgeons preferred the GS solution due to a smaller bone protrusion with respect to the corresponding OA solution. In the case shown in Figure 12B, the GS solution has a better plate positioning, albeit a slightly worse reduction alignment compared to the OA solution. In 8 out of the 36 cases, the GS solution corresponded to a closing wedge osteotomy (C1, C5-C7, C10, C18, C20, and C21) because planners failed to manually find a clinically acceptable single-cut osteotomy. For all of these 8 cases, the optimization framework was rerun and the algorithm converged towards a single-cut solution by applying the corresponding constraints. A single-cut osteotomy is the most preferred clinical implementation of a corrective osteotomy. Figure 13 shows one of the OA single-cut solutions in comparison to the closing- wedge approach of the GS (Figure 13A). In the OA solution (Figure 13B), not only the osteotomy type was improved, but also the position of the fixation plate. Figure 11: Example cases for OA solutions that outperformed GS solutions. The proximal fragment of the pathological bone is shown in white and aligned to the reconstruction target, shown in green. The reduced distal fragments of the OA and GS solutions are shown in magenta and purple, respectively. (A) Case 21 depicts a radius shaft osteotomy where the OA solution was able to solve undesired plate-bone intersections. (B) Case 30 shows a distal radius osteotomy. The OA solution was able to find a better positioning of the fixation plate avoiding bone intersections and improving the bone-plate contact. Figure 12: Example cases where GS solutions were assessed to be better than OA solutions. The proximal fragment of the pathological bone is shown in white and aligned to the reconstruction target, shown in green. The reduced distal fragments of the OA and GS solutions are shown in magenta and purple, respectively. (A) Case 07 depicts a radius shaft osteotomy where the bone protrusion was assessed to be better in the GS solution. The bone protrusion in the OA solution was larger in comparison to the GS, albeit a slightly better reduction alignment (B) Case 26 shows a distal radius osteotomy where the position of the plate was assessed to be better in the GS solution. In the OA solution, the plate-bone distance was optimized for the proximal fragment, resulting in a small bone intersection on the distal fragment. Figure 13: Example of case 20, where the algorithm was able to find a single-cut osteotomy solution. The proximal fragment of the pathological bone is shown in white and the reconstruction target is shown in green. (A) Closing-wedge solution generated by the state-of-the-art planning (GS); realigned distal fragment is shown in magenta and the wedge to be removed to complete the osteotomy is shown in red. (B) Single-cut solution generated by our approach (OA); realignment of the distal fragment is shown in purple. Position of the fixation plate was also improved with respect to the GS solution. 3.3. Quantitative Evaluation Six different error measurements were evaluated for both, the GS and OA solutions, and for each of the 36 cases. Their average values, standard deviations, and ranges are reported in Table 5. (a) Reduction error, measured by Eq. 6, between the transformed distal fragment and the reconstruction target. A reduction error <= 0.5 mm is desired clinically. (b) Average distance between the fixation plate and the proximal bone fragment measured by Eq. 9. (c) Average distance between the fixation plate and the distal bone fragment measured by Eq. 10. (d) Average screw purchase of the proximal screws of the fixation plate, measured by Eq. 14. A larger screw purchase represents a better positioning of the screw. (e) Average screw purchase of the distal screws of the fixation plate, measured by Eq. 14. A larger screw purchase represents a better positioning of the screw. (f) Minimum distance between the screws of the fixation plate and the osteotomy plane measured by 𝑚𝑚𝑒𝑒𝑛𝑛 � 𝑛𝑛� ⃗ . (𝐶𝐶 − 𝑙𝑙𝑝𝑝 )� � � 𝑛𝑛� ⃗ � . A distance ≥2 mm is preferred. 𝑃𝑃 𝑝𝑝𝑙𝑙 𝑚𝑚 𝑝𝑝𝑙𝑙 𝑖𝑖 Table 5: Mean, standard deviation and range for the 6 error measurements used in the technical evaluation of the optimization solutions compared to the gold standard. Gold Standard (GS) Optimization Algorithm (OA) Measurement Mean ± SD Range Mean ± SD Range Alignment Error (mm) 1.47 1.89 0.08 , 9.56 0.94 0.88 0.07 , 3.55 Distance Plate-Bone Proximal (mm) 0.83 0.67 -1.04 , 2.64 0.57 0.44 -0.80 , 1.59 Distance Plate-Bone Distal (mm) 0.88 0.66 -1.05 , 2.03 0.64 0.44 -0.88 , 1.62 Screw Purchase Proximal (unit) 1.4 1.02 -0.86 , 2.90 1.98 0.74 -0.22 , 3.34 Screw Purchase Distal (unit) 1.15 0.89 -0.22 , 3.26 1.32 0.73 -0.06 , 2.69 Minimum Distance 4.91 3.05 0.45 , 9.59 5.75 2.53 0.80 , 12.34 Screw-Osteotomy Plane (mm) We have also performed an evaluation with respect to each planning objective by assessing the error measures for both GS and OA solutions. Figure 14 shows a set of box plots expressed as the relative improvement of the measures that were achieved by the OA solution. The blue baseline represents the GS value. In average, the alignment error was improved by 20% with respect to the GS among all the evaluated cases. The distance between the plate and the proximal and distal bone fragments has been improved by 26% and 31%, respectively. The average screw purchase for the proximal and distal screws was improved by 106% and 107%, respectively. Finally, an improvement of 93% in the average minimum distance between the fixation screws and the osteotomy plane was achieved by the OA solutions. Figure 14: Results of quantitative evaluation. Box plot (3/2 interquartile range whiskers) of the relative improvement of OA solutions in relation to the GS (blue vertical bold line) for each of the 36 cases. The black bold line indicates the median and the red star indicates the average improvement for each error measurement. The black crosses indicate the outliers. 3.4. Runtime The presented approach was implemented in Matlab R2017b (MATLAB, 2017). The average runtime of the algorithm comprising pre-processing, automatic diagnosis, verification and initialization of constraints, two optimization stages and generation of output solution was 85.38 min (±33.15 min). 4. Discussion In the last decade, the use of computer-assisted 3D preoperative planning for orthopedic surgeries has increased significantly due to its higher precision and to the capability of treating more complex pathologies (Murase et al., 2008; Nagy et al., 2008). However, the great effort required for generating preoperative planning solutions with current state-of-the-art approaches poses a bottleneck in the treatment of corrective osteotomies. In this work, we have presented a computer framework for the fully automatic calculation of preoperative planning solutions of osteotomies. Our approach augments the skills of the surgeons by providing improved solutions with a reliable computational tool. To the best of our knowledge, this automated preoperative planning framework represents the first to have been developed to generate ready-to-use preoperative planning solutions for corrective osteotomies of the forearm In orthopedic surgery, good surgical outcome is tightly associated with careful preoperative planning and precise surgical execution (Miyake et al., 2012b; Vlachopoulos et al., 2015; Vroemen et al., 2012). In order to demonstrate the clinical validity of our framework, we conducted a clinical study which confirmed the capabilities of the algorithm for generating ready-to-use preoperative solutions. The study showed that solutions generated by our algorithm were preferred in 55% of the time. In the remaining 45%, the quantitative evaluation confirmed still a high quality of the OA solution, with an error to the GS solutions of less than 1 mm for the reduction alignment, and less than 1.5 mm for the bone-plate distance. One clear advantage of our framework is the ability of generating solutions in an automated fashion, deducting human calculation and planning times from the clinical setting, in contrast to state-of-the-art methods, which are associated with high effort and clinical costs (Fürnstahl et al., 2016). The calculation of 3D preoperative planning solutions in an automatic way would represent a time reduction of almost 35% of the human workload per case, decreasing therefore significantly the related costs. Our institution performs more than 250 corrective osteotomies each year using a commercial service, in which the 3D preoperative planning is performed in a manual fashion. If we consider average diagnosis and planning costs of USD 600 per case, our approach would save USD 150’000 of yearly clinical costs in our institution alone, while improving also the quality of planning and patient treatment. Another advantage is the capability of calculating solutions with different trade-offs among the optimization objectives. Our algorithm can drastically improve one of the objectives by slightly decreasing the quality of other objectives. Such evaluations are not possible to obtain in a systematic way using the manual state-of-the-art approach (Bauer et al., 2017b; Vlachopoulos et al., 2015). Moreover, our method allows calculating single-cut osteotomies for arbitrary bone deformities including angular and translational deformities. In the presented case series, 8 single-cut OA solutions could be calculated where the real surgery was implemented with a closing wedge osteotomy. In 75% of these cases, readers validated the generated OA single-cut osteotomy solutions as a better preoperative planning than the original closing-wedge osteotomy proposed by a surgeon. Concerning the surgical execution, the generation of complex solutions for preoperative planning comes also at the cost of an increase in the complexity of the intraoperative navigation. At our institution, we have developed a state-of-the-art technique for intraoperative guidance using PSI, which are 3D-printed individually for each case. These PSI allow a high level of precision for executing the preoperative plan in the surgery (Fürnstahl et al., 2016; Omori et al., 2014; Rosseels et al., 2018). An example of the generated PSI for intraoperative navigation of a distal radius osteotomy is shown in Figure 15. The accuracy of PSI navigation for corrective osteotomies has been previously evaluated, showing an improve in the postoperative outcome (Asada et al., 2014; Bauer et al., 2017a; Byrne et al., 2017; Caiti et al., 2018; Farshad et al., 2017; Gouin et al., 2014; Kunz et al., 2013; Michielsen et al., 2018; Omori et al., 2014; Roner et al., 2018; Rosseels et al., 2018; Schweizer et al., 2016; Vlachopoulos et al., 2015) Figure 15: Example of translation of preoperative plan to the intraoperative situation by means of 3D-printed PSI. (A) Preoperative plan of a distal radius osteotomy. (B) PSI design for the osteotomy cut. (C) PSI design for the implant navigation and bone reduction. (D) Intraoperative navigation of the osteotomy cut. (E) Intraoperative navigation of the bone reduction and implant positioning. The black arrows indicate the PSI design that corresponds to each intraoperative situation. The presented framework has some technical and clinical limitations. In the first place, the number of cases included in this study was not enough to provide statistical relevance to the quantitative and qualitative analysis performed. However, the observed trend shows a promising potential of the capabilities of the developed framework. The quality of the solutions generated by the optimization algorithm was very similar for every case to the ones of the Gold Standard. The decision towards one or the other was determined by very specific factors such as a better angle for the osteotomy cut or a better position for the fixation plate. Nevertheless, this is also the case when treating patients with the Gold Standard method as the surgeon has to decide which solution will be the one transferred to the intraoperative situation by means of PSI (Roner et al., 2018; Rosseels et al., 2018; Schweizer et al., 2013; Victor and Premanathan, 2013; Vlachopoulos et al., 2015; Wong et al., 2017) In order to study the accuracy of surgical outcomes using automatic optimization methods, a prospective study with a statistical significant number of cases needs to be performed including a control group. Nevertheless, before such a clinical study can be approved and performed, the performance and safety of the algorithm needed to be validated. We have successfully achieved that in this paper through two different steps: (1) a quantitative evaluation showing that the algorithm is able to generate the expected solutions, optimizing all objectives, and (2) qualitative evaluation by a survey among surgeons and engineers showing that the algorithm produces solutions that the operating surgeon would execute in the practice. Secondly, the optimization is not able to take into account clinical decisions made by the surgeons based on patient history, for example, an unusual reduction alignment due to soft tissue pathologies (Vlachopoulos et al., 2015), or a different fixation technique due to osteoporotic bones. Also, the presence of a deformity in the surrounding bones has a significant impact in the final reduction. In our algorithm, this influence was taken into account through the reconstruction target, however in future works, the influence of neighboring bones should be also taken into account as a part of the optimization framework. Furthermore, surgeons have intraoperative techniques to improve poorly fitting plates such as manual plate-bending or bone reshaping. Such actions cannot be foreseen by the algorithm. This might be also a reason why the GS solution was preferred in some cases even with a poor fitting of the fixation plate. Example of that can be seen for cases 7, 17, and 20, where the GS solution was preferred by 83% of the readers. In the future, patient-specific osteosynthesis and osteotomy implants such as the one described by Caiti et al. (2018) and Omori et al. (2014) could be integrated into the optimization framework. Our approach is an important first step to turn automatic osteotomy planning into clinical practice. Due to the modularity of the algorithm, the approach could likely be extended to handle further anatomies (e.g., hip, humerus, clavicle, wrist bones) (Bugeau and Pérez, 2008; Chahla et al., 2016; Reising et al., 2013; Tschannen et al., 2016; Vlachopoulos et al., 2018). The positioning of the fixation plate and the calculation of the screw purchase could directly be adapted without significant efforts, but the calculation of the osteotomy plane and bone reduction would require the adaptation of the fitness functions and constraints. Our algorithm does not include yet the possibility to support simultaneous osteotomies on multiple locations (Belei et al., 2007; Fürnstahl et al., 2016; Schkommodau et al., 2005), and it does not include either the capability to plan intra-articular osteotomies (Schweizer et al., 2013). However, we consider that the strategy could be adapted by including extra stages to the optimization pipeline, at the cost of increasing the calculation times due to the linear structure of the optimization stages. One possibility to avoid the increase of calculation times is to implement a parallel design for the optimization stages, also allowing the capability for a back- feedback mechanism of the entire optimization. The presented approach is 18% faster in comparison to our previous work (Carrillo et al., 2017). However, the computational effort is still very high. One of the reasons for long calculation times is the implementation of the algorithm in a non-compiled programming language. Another possible source for the long calculation time is the lack of parallelization. We expect that the implementation of the entire pipeline in a compiled language and the use of graphical processing units (GPUs) for a customized parallelization of the approach will speed up the optimization framework to interactive rates. 5. Conclusions and Further work The presented framework is able to generate clinically feasible preoperative planning solutions in an automatic fashion. The key idea of the approach is the use of a multi-staged, multi-objective optimization pipeline for the concurrent optimization of the planning objectives. We demonstrated that the clinical and technical quality of the solutions generated by the algorithm can be of the same quality as the solutions created by experienced surgeons. Moreover, in more than 50% of the cases, our algorithm was able to generate solutions that are hardly possible to be found within a reasonable time by the trial and error method of the Gold Standard. This shows the potential of the method to improve the quality of the preoperative planning by considering a larger range of feasible preoperative solutions. The algorithm does not intend to replace the expertise of surgeons in the process; our aim is to reduce unnecessary human work load. We seek to provide surgeons with automated tools for patient diagnosis and treatment, which can substantially decrease associated clinical times and costs, and generate better clinical outputs. Future works will focus on more flexible and less-invasive techniques for intraoperative guidance than patient-specific instrument (PSI) approaches such as augmented reality techniques in combination with learning-based methods for more comprehensive surgical guidance. Future solutions could also include a completely adaptive framework, were the optimal preoperative planning solution is adapted in real time according to the intraoperative situation. Conflict of interest None of the authors have any financial or personal relationships with other people or organizations that could bias their work. Acknowledgements This work has been funded through a Balgrist Foundation grant and a highly specialized medicine grant (HSM2) of the Canton of Zurich. We would also like to acknowledge the support and valuable input from the surgeons and planning engineers at the department of orthopaedics of Balgrist University Hospital. References Andress, S., Johnson, A., Unberath, M., Winkler, A.F., Yu, K., Fotouhi, J., Weidert, S., Osgood, G., Navab, N., 2018. On- the-fly augmented reality for orthopedic surgery using a multimodal fiducial. J Med Imaging (Bellingham) 5, 021209. Asada, S., Mori, S., Matsushita, T., Nakagawa, K., Tsukamoto, I., Akagi, M., 2014. Comparison of MRI- and CT-based patient-specific guides for total knee arthroplasty. The Knee 21, 1238-1243. Athwal, G.S., Ellis, R.E., Small, C.F., Pichora, D.R., 2003. Computer-assisted distal radius osteotomy1. The Journal of Hand Surgery 28, 951-958. Bauer, A.S., Storelli, D.A.R., Sibbel, S.E., Mccarroll, H.R., Lattanza, L.L., 2017a. Preoperative computer simulation and patient-specific guides are safe and effective to correct forearm deformity in children. Journal of Pediatric Orthopaedics 37, 504-510. Bauer, D.E., Zimmermann, S., Aichmair, A., Hingsammer, A., Schweizer, A., Nagy, L., Fürnstahl, P., 2017b. Conventional Versus Computer-Assisted Corrective Osteotomy of the Forearm: a Retrospective Analysis of 56 Consecutive Cases. The Journal of Hand Surgery 42, 447-455. Belei, P., Schkommodau, E., Frenkel, A., Mumme, T., Radermacher, K., 2007. Computer-assisted single- or double-cut oblique osteotomies for the correction of lower limb deformities. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 221, 787-800. Besl, P.J., McKay, N.D., 1992. Method for registration of 3-D shapes, Sensor Fusion IV: Control Paradigms and Data Structures. International Society for Optics and Photonics, pp. 586-607. Bilic, R., Zdravkovic, V., Boljevic, Z., 1994. Osteotomy for deformity of the radius. Computer-assisted three-dimensional modelling. Journal of Bone & Joint Surgery, British Volume 76-B, 150-154. Blecha, L.D., Zambelli, P.Y., Ramaniraka, N.A., Bourban, P.E., Månson, J.A., Pioletti, D.P., 2005. How plate positioning impacts the biomechanics of the open wedge tibial osteotomy; A finite element analysis. Computer Methods in Biomechanics and Biomedical Engineering 8, 307-313. Bugeau, A., Pérez, P., 2008. Track and Cut: Simultaneous Tracking and Segmentation of Multiple Objects with Graph Cuts. EURASIP Journal on Image and Video Processing 2008, 317278. Byrne, A.-M., Impelmans, B., Bertrand, V., Van Haver, A., Verstreken, F., 2017. Corrective Osteotomy for Malunited Diaphyseal Forearm Fractures Using Preoperative 3-Dimensional Planning and Patient-Specific Surgical Guides and Implants. Journal of Hand Surgery 42, 836.e831-836.e812. Caiti, G., Dobbe, J.G.G., Loenen, A.C.Y., Beerens, M., Strackee, S.D., Strijkers, G.J., Streekstra, G.J., 2018. Implementation of a semiautomatic method to design patient-specific instruments for corrective osteotomy of the radius. International Journal of Computer Assisted Radiology and Surgery. Carrillo, F., Vlachopoulos, L., Schweizer, A., Nagy, L., Snedeker, J., Fürnstahl, P., 2017. A Time Saver: Optimization Approach for the Fully Automatic 3D Planning of Forearm Osteotomies, Medical Image Computing and Computer-Assisted Intervention − MICCAI 2017. Springer International Publishing, Quebec, Canada, pp. 488-496. Chahla, J., Dean, C.S., Mitchell, J.J., Moatshe, G., Serra Cruz, R., LaPrade, R.F., 2016. Medial Opening Wedge Proximal Tibial Osteotomy. Arthroscopy Techniques 5, e919-e928. Criminisi, A., Shotton, J., 2013. Decision forests for computer vision and medical image analysis. Springer Science & Business Media. Deb, K., Agrawal, S., Pratap, A., Meyarivan, T., 2000. A fast elitist non-dominated sorting genetic algorithm for multi- objective optimization: NSGA-II. Lecture notes in computer science 1917, 849-858. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. Evolutionary Computation, IEEE Transactions on 6, 182-197. Dobbe, J.G.G., Strackee, S.D., Schreurs, A.W., Jonges, R., Carelsen, B., Vroemen, J.C., Grimbergen, C.A., Streekstra, G.J., 2011. Computer-Assisted Planning and Navigation for Corrective Distal Radius Osteotomy, Based on Pre- and Intraoperative Imaging. Biomedical Engineering, IEEE Transactions on 58, 182-190. Dobbe, J.G.G., Strackee, S.D., Streekstra, G.J., 2017. Minimizing the translation error in the application of an oblique single- cut rotation osteotomy: Where to cut? IEEE Transactions on Biomedical Engineering PP, 1-1. Dobbe, J.G.G., Vroemen, J.C., Strackee, S.D., Streekstra, G.J., 2014. Patient-specific distal radius locking plate for fixation and accurate 3D positioning in corrective osteotomy. Strategies in Trauma and Limb Reconstruction 9, 179-183. Dougherty, P.J., Kim, D.-G., Meisterling, S., Wybo, C., Yeni, Y., 2008. Biomechanical Comparison of Bicortical Versus Unicortical Screw Placement of Proximal Tibia Locking Plates: A Cadaveric Model. Journal of Orthopaedic Trauma 22, 399- Eck, M., Hoschek, J., Weber, U., 1990. Three-dimensional determination of an oblique osteotomy in the hip by mathematical optimization fulfilling some anatomical demands. Journal of biomechanics 23, 1061-1067. Esfandiari, H., Newell, R., Anglin, C., Street, J., Hodgson, A.J., 2018. A deep learning framework for segmentation and pose estimation of pedicle screw implants based on C-arm fluoroscopy. International Journal of Computer Assisted Radiology and Surgery 13, 1269-1282. Fang, Q., Boas, D.A., 2009. Tetrahedral mesh generation from volumetric binary and grayscale images, Biomedical Imaging: From Nano to Macro, 2009. ISBI'09. IEEE International Symposium on. Ieee, pp. 1142-1145. Farshad, M., Betz, M., Farshad-Amacker, N.A., Moser, M., 2017. Accuracy of patient-specific template-guided vs. free-hand fluoroscopically controlled pedicle screw placement in the thoracic and lumbar spine: a randomized cadaveric study. European Spine Journal 26, 738-749. Fernandez, D.L., 1982. Correction of post-traumatic wrist deformity in adults by osteotomy, bone-grafting, and internal fixation. J Bone Joint Surg Am 64, 1164-1178. Fürnstahl, P., 2010. Computer-assisted planning for orthopedic surgery. Diss., Eidgenössische Technische Hochschule ETH Zürich, Nr. 19102, 2010. Fürnstahl, P., Fuchs, T., Schweizer, A., Nagy, L., Székely, G., Harders, M., 2008. Automatic and robust forearm segmentation using graph cuts, Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on. IEEE, pp. 77-80. Fürnstahl, P., Schweizer, A., Graf, M., Vlachopoulos, L., Fucentese, S., Wirth, S., Nagy, L., Szekely, G., Goksel, O., 2016. Surgical Treatment of Long-Bone Deformities: 3D Preoperative Planning and Patient-Specific Instrumentation, In: Zheng, G., Li, S. (Eds.), Computational Radiology for Orthopaedic Interventions. Springer International Publishing, Cham, pp. 123- Gosse, F., Brack, C., Götte, H., Roth, M., Rühmann, O., Schweikard, A., Vahldiek, M., 1997. Roboterunterstützung in der KnieendoprothetikRobotic-assisted total knee replacement. Der Orthopäde 26, 258-266. Gouin, F., Paul, L., Odri, G.A., Cartiaux, O., 2014. Computer-Assisted Planning and Patient-Specific Instruments for Bone Tumor Resection within the Pelvis: A Series of 11 Patients. Sarcoma 2014, 1-9. Jolliffe, I., 2011. Principal component analysis, International encyclopedia of statistical science. Springer, pp. 1094-1096. Kawakami, H., Sugano, N., Nagaoka, T., Hagio, K., Yonenobu, K., Yoshikawa, H., Ochi, T., Hattori, A., Suzuki, N., 2002. 3D analysis of the alignment of the lower extremity in high tibial osteotomy, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, pp. 261-267. Klausmeyer, M.A., Mudgal, C., 2014. Exposure of the Forearm and Distal Radius. Hand Clinics 30, 427-433. Kubiak, E.N., Fulkerson, E., Strauss, E., Egol, K.A., 2006. The evolution of locked plates. JBJS 88, 189-200. Kunz, M., Ma, B., Rudan, J.F., Ellis, R.E., Pichora, D.R., 2013. Image-Guided Distal Radius Osteotomy Using Patient- Specific Instrument Guides. The Journal of Hand Surgery 38, 1618-1624. Lorensen, W.E., Cline, H.E., 1987. Marching cubes: A high resolution 3D surface construction algorithm, ACM siggraph computer graphics. ACM, pp. 163-169. Lüthi, M., Gerig, T., Jud, C., Vetter, T., 2018. Gaussian Process Morphable Models. IEEE Transactions on Pattern Analysis and Machine Intelligence 40, 1860-1873. MATLAB, 2017. 9.3.0.713579 (R2017b). The MathWorks Inc., Natick, Massachusetts. Mauler, F., Langguth, C., Schweizer, A., Vlachopoulos, L., Gass, T., Lüthi, M., Fürnstahl, P., 2017. Prediction of Normal Bone Anatomy for the Planning of Corrective Osteotomies of Malunited Forearm Bones Using a Three-Dimensional Statistical Shape Model. Journal of orthopaedic research: official publication of the Orthopaedic Research Society 35, 2630- Menetrey, J., Paul, M., 2004. Möglichkeiten der computergestützten Navigation bei kniegelenknahen Osteotomien. Der Orthopäde 33, 224-228. Merle d’Aubigné, R., Descamps, L., 1952. L’ostéotomie plane oblique dans la correction des deformations des membres. Bull Mem Arch Chirurgie 8, 271-276. Michielsen, M., Van Haver, A., Bertrand, V., Vanhees, M., Verstreken, F., 2018. Corrective osteotomy of distal radius malunions using three-dimensional computer simulation and patient-specific guides to achieve anatomic reduction. European journal of orthopaedic surgery & traumatology : orthopedie traumatologie 28, 1531-1535. Miettinen, K., 1999. Nonlinear Multiobjective Optimization. Springer US. Miyake, J., Moritomo, H., Kataoka, T., Murase, T., Sugamoto, K., 2012a. In vivo three-dimensional motion analysis of chronic radial head dislocations. Clinical Orthopaedics and Related Research® 470, 2746-2755. Miyake, J., Murase, T., Moritomo, H., Sugamoto, K., Yoshikawa, H., 2011. Distal Radius Osteotomy with Volar Locking Plates Based on Computer Simulation. Clinical Orthopaedics and Related Research® 469, 1766-1773. Miyake, J., Murase, T., Oka, K., Moritomo, H., Sugamoto, K., Yoshikawa, H., 2012b. Three-Dimensional Corrective Osteotomy for Malunited Diaphyseal Forearm Fractures Using Custom-Made Surgical Guides Based on Computer Simulation. JBJS Essential Surgical Techniques 2, e24. Mount, D.M., Arya, S., 1998. ANN: library for approximate nearest neighbour searching. Murase, T., Oka, K., Moritomo, H., Goto, A., Yoshikawa, H., Sugamoto, K., 2008. Three-Dimensional Corrective Osteotomy of Malunited Fractures of the Upper Extremity with Use of a Computer Simulation System. The Journal of Bone & Joint Surgery 90, 2375-2389. Nagy, L., Jankauskas, L., Dumont, C.E., 2008. Correction of forearm malunion guided by the preoperative complaint. Clinical orthopaedics and related research 466, 1419-1428. Omori, S., Murase, T., Kataoka, T., Kawanishi, Y., Oura, K., Miyake, J., Tanaka, H., Yoshikawa, H., 2014. Three‐ dimensional corrective osteotomy using a patient‐ specific osteotomy guide and bone plate based on a computer simulation system: accuracy analysis in a cadaver study. The International Journal of Medical Robotics and Computer Assisted Surgery 10, 196-202. Pineda, J., 1988. A parallel algorithm for polygon rasterization, ACM SIGGRAPH Computer Graphics. ACM, pp. 17-20. Reising, K., Strohm, P.C., Hauschild, O., Schmal, H., Khattab, M., Südkamp, N.P., Niemeyer, P., 2013. Computer-assisted navigation for the intraoperative assessment of lower limb alignment in high tibial osteotomy can avoid outliers compared with the conventional technique. Knee Surgery, Sports Traumatology, Arthroscopy 21, 181-188. Reitman, C., Nguyen, L., Fogel, G., 2004. Biomechanical evaluation of relationship of screw pullout strength, insertional torque, and bone mineral density in the cervical spine. Journal of spinal disorders & techniques 17, 306-311. Roner, S., Carrillo, F., Vlachopoulos, L., Schweizer, A., Nagy, L., Fuernstahl, P., 2018. Improving accuracy of opening- wedge osteotomies of distal radius using a patient-specific ramp-guide technique. BMC musculoskeletal disorders 19, 374. Roner, S., Vlachopoulos, L., Nagy, L., Schweizer, A., Fürnstahl, P., 2017. Accuracy and Early Clinical Outcome of 3- Dimensional Planned and Guided Single-Cut Osteotomies of Malunited Forearm Bones. The Journal of Hand Surgery 42, 1031.e1031-1031.e1038. Roscoe, L., 1988. Stereolithography interface specification. America-3D Systems Inc 27. Rosseels, W., Herteleer, M., Sermon, A., Nijs, S., Hoekstra, H., 2018. Corrective osteotomies using patient-specific 3D- printed guides: a critical appraisal. European Journal of Trauma and Emergency Surgery. Sangeorzan, B.P., Judd, R.P., Sangeorzan, B.J., 1989. Mathematical analysis of single-cut osteotomy for complex long bone deformity. Journal of Biomechanics 22, 1271-1278. Schenk, P., Vlachopoulos, L., Hingsammer, A., Fucentese, S.F., Fürnstahl, P., 2016. Is the contralateral tibia a reliable template for reconstruction: a three-dimensional anatomy cadaveric study. Knee Surg Sports Traumatol Arthrosc. Scheufler, K.-M., Franke, J., Eckardt, A., Dohmen, H., 2011. Accuracy of Image-Guided Pedicle Screw Placement Using Intraoperative Computed Tomography-Based Navigation With Automated Referencing, Part I: Cervicothoracic Spine. Neurosurgery 69, 782-795. Schkommodau, E., Frenkel, A., Belei, P., Recknagel, B., Wirtz, D.-C., Radermacher, K., 2005. Computer-assisted optimization of correction osteotomies on lower extremities. Computer Aided Surgery 10, 345-350. Schneider, P., Eberly, D.H., 2002. Geometric tools for computer graphics. Elsevier. Schweizer, A., Fürnstahl, P., Harders, M., Székely, G., Nagy, L., 2010. Complex radius shaft malunion: osteotomy with computer-assisted planning. Hand 5, 171-178. Schweizer, A., Fürnstahl, P., Nagy, L., 2013. Three-dimensional correction of distal radius intra-articular malunions using patient-specific drill guides. The Journal of hand surgery 38, 2339-2347. Schweizer, A., Mauler, F., Vlachopoulos, L., Nagy, L., Fürnstahl, P., 2016. Computer-assisted 3-dimensional reconstructions of scaphoid fractures and nonunions with and without the use of patient-specific guides: early clinical outcomes and postoperative assessments of reconstruction accuracy. The Journal of hand surgery 41, 59-69. Sepasian, N., Van de Giessen, M., Dobbe, I., Streekstra, G., 2014. Bone reposition planning for corrective surgery using statistical shape model: assessment of differential geometrical features, Bayesian and grAphical Models for Biomedical Imaging. Springer, pp. 49-60. Spahn, G., 2004. Complications in high tibial (medial opening wedge) osteotomy. Archives of Orthopaedic and Trauma Surgery 124, 649-653. Subburaj, K., Ravi, B., Agarwal, M., 2010. Computer-aided methods for assessing lower limb deformities in orthopaedic surgery planning. Computerized Medical Imaging and Graphics 34, 277-288. Tschannen, M., Vlachopoulos, L., Gerber, C., Székely, G., Fürnstahl, P., 2016. Regression forest-based automatic estimation of the articular margin plane for shoulder prosthesis planning. Medical Image Analysis 31, 88-97. Vatti, B.R., 1992. A generic solution to polygon clipping. Communications of the ACM 35, 56-63. Victor, J., Premanathan, A., 2013. Virtual 3D planning and patient specific surgical guides for osteotomies around the knee: a feasibility and proof-of-concept study. Bone Joint J 95-B, 153-158. Vlachopoulos, L., Carrillo, F., Gerber, C., Székely, G., Fürnstahl, P., 2017. A Novel Registration-Based Approach for 3D Assessment of Posttraumatic Distal Humeral Deformities. JBJS 99, e127. Vlachopoulos, L., Schweizer, A., Graf, M., Nagy, L., Fürnstahl, P., 2015. Three-dimensional postoperative accuracy of extra- articular forearm osteotomies using CT-scan based patient-specific surgical guides. BMC musculoskeletal disorders 16, 1. Vlachopoulos, L., Székely, G., Gerber, C., Fürnstahl, P., 2018. A scale-space curvature matching algorithm for the reconstruction of complex proximal humeral fractures. Medical Image Analysis 43, 142-156. Vroemen, J.C., Dobbe, J.G.G., Jonges, R., Strackee, S.D., Streekstra, G.J., 2012. Three-Dimensional Assessment of Bilateral Symmetry of the Radius and Ulna for Planning Corrective Surgeries. The Journal of Hand Surgery 37, 982-988. Wall, L.B., Brodt, M.D., Silva, M.J., Boyer, M.I., Calfee, R.P., 2012. The Effects of Screw Length on Stability of Simulated Osteoporotic Distal Radius Fractures Fixed With Volar Locking Plates. The Journal of Hand Surgery 37, 446-453. Weninger, P., Dall'Ara, E., Leixnering, M., Pezzei, C., Hertz, H., Drobetz, H., Redl, H., Zysset, P., 2010. Volar fixed-angle plating of extra-articular distal radius fractures—a biomechanical analysis comparing threaded screws and smooth pegs. Journal of Trauma and Acute Care Surgery 69, E46-E55. Wong, K., Paul, L., Gerbers, J., 2017. PATIENT SPECIFIC INSTRUMENTS CAN ACHIEVE A BETTER SURGICAL ACCURACY THAN NAVIGATION ASSISTANCE IN JOINT-PRESERVING SURGERY OF THE KNEE JOINT: A CADAVERIC COMPARATIVE STUDY. CAOS 1, 104-108. Zdravkovic, V., Bilic, R., 1990. Computer-assisted preoperative planning (CAPP) in orthopaedic surgery. Computer Methods and Programs in Biomedicine 32, 141-146. Vitae: Fabio Carrillo received his master degree in electronic engineering from Universidad Simon Bolivar in Caracas, Venezuela, with focus on automation. In 2015, he received his MSc degree in Robotics, Systems and Control from ETH Zürich, with a deep research interest in biomedical engineering. He joined the Computer Assisted Research and Development (CARD) Group at the University Hospital Balgrist in Zurich in September 2015, where he performs his doctoral studies, together with the Laboratory for Orthopaedic Biomechanics of the Department of Health Sciences and Technology at the ETH Zürich. Simon Roner is an orthopaedic surgeon resident working as a research associate in computer-assisted surgery. He graduated from the medical school at Zurich University in 2012. After completion of Swiss board exams in orthopaedic surgery, he will pursue a specialization in hand and upper limb surgery. Marco von Atzigen received his Master’s degree in Mechanical Engineering from ETH Zurich with focus on Robotics, Systems and Control in July 2017. He joined CARD in October 2017 as a software developer working on machine learning and augmented reality. Since December 2018, he is pursuing his doctoral studies in deep learning for medical imaging and scene understanding between CARD, the spine surgery team at the University Hospital Balgrist Zürich and the Laboratory for Orthopaedic Biomechanics of the Department of Health Sciences and Technology at the ETH Zürich. Andreas Schweizer graduated in Medicine from the University of Zurich. He received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 2003 and of Hand Surgery in 2005. He completed his MD thesis at Institute of Anatomy, University of Bern and his state doctorate at the University of Zurich. He is currently working as a consultant hand surgeon at the Balgrist University Hospital in Zurich, Switzerland and is specialized in 3D planning of orthopedic surgeries of upper extremity. Ladislav Nagy graduated in Medicine from the University of Zurich in 1982. He received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 1990 and of Hand Surgery in 1992. He is titular professor of the faculty of medicine of the University of Zurich since 2011. He is the current chief of the hand surgery of the university hospital Balgrist and work as a consultant specialist in 3D planning of orthopedic surgeries of the hand. Lazaros Vlachopoulos graduated in Medicine from the RWTH Aachen in Germany in 2004 and completed his MD thesis at the Institute of Human Genetics, RWTH Aachen. Afterwards, he completed his residency in orthopaedic surgery in Germany and Switzerland and received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 2012. In 2017 he received his PhD in medical image analysis from the ETH Zurich, Switzerland. He is currently working as a consultant orthopaedic surgeon at the Balgrist University Hospital in Zurich, Switzerland, specialized in computer-assisted orthopaedic surgery. Jess Snedeker is Associate Professor of Orthopaedic Biomechanics, holding a professorial chair at both the ETH Zurich (Department of Health Sciences and Technology) and the Medical Faculty of the University of Zurich (Department of Orthopaedics). He heads the division of experimental research at the University Hospital Balgrist, and also serves as the chief scientific officer of the Balgrist Campus, designated in 2017 by the Swiss Secretariat for Education, Research, and Innovation as a “Research Infrastructure of National Relevance”. Philipp Fürnstahl received the MSc degree in technical mathematics and information procession from the Technical University of Graz, Austria, in 2005. In 2010, he received the PhD degree in medical image analysis from the ETH Zurich, Switzerland, for his research in the field of computer-assisted preoperative surgery planning. Philipp Fürnstahl is currently head of the Computer Assisted Research and Development (CARD) Group at the University Hospital Balgrist in Zurich, Switzerland, where he is involved in the development of patient-specific solutions for orthopaedics surgeries. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

An Automatic Genetic Algorithm Framework for the Optimization of Three-dimensional Surgical Plans of Forearm Corrective Osteotomies

Loading next page...
 
/lp/arxiv-cornell-university/an-automatic-genetic-algorithm-framework-for-the-optimization-of-three-ychUKIEfpq

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
1361-8415
eISSN
ARCH-3348
DOI
10.1016/j.media.2019.101598
Publisher site
See Article on Publisher Site

Abstract

Three-dimensional (3D) computer-assisted corrective osteotomy has become the state-of-the-art for surgical treatment of complex bone deformities. Despite available technologies, the automatic generation of clinically acceptable, ready-to-use preoperative planning solutions is currently not possible for such pathologies. Multiple contradicting and mutually dependent objectives have to be considered, as well as clinical and technical constraints, which generally require iterative manual adjustments. This leads to unnecessary surgeon efforts and unbearable clinical costs, hindering also the quality of patient treatment due to the reduced number of solutions that can be investigated in a clinically acceptable timeframe. In this paper, we propose an optimization framework for the generation of ready-to-use preoperative planning solutions in a fully automatic fashion. An automatic diagnostic assessment using patient-specific 3D models is performed for 3D malunion quantification and definition of the optimization parameters’ range. Afterward, clinical objectives are translated into the optimization module, and controlled through tailored fitness functions based on a weighted and multi-staged optimization approach. The optimization is based on a genetic algorithm capable of solving multi-objective optimization problems with non-linear constraints. The framework outputs a complete preoperative planning solution including position and orientation of the osteotomy plane, transformation to achieve the bone reduction, and position and orientation of the fixation plate and screws. A qualitative validation was performed on 36 consecutive cases of radius osteotomy where solutions generated by the optimization algorithm (OA) were compared against the gold standard solutions generated by experienced surgeons (Gold Standard; GS). Solutions were blinded and presented to 6 readers (4 surgeons, 2 planning engineers), who voted OA solutions to be better in 55% of the time. The quantitative evaluation was based on different error measurements, showing average improvements with respect to the GS from 20% for the reduction alignment and up to 106% for the position of the fixation screws. Notably, our algorithm was able to generate feasible clinical solutions which were not possible to obtain with the current state-of-the-art method. Keywords 3D Surgical Planning · Automatic · Forearm · Corrective Osteotomy · Multi-objective Optimization Research Highlights • Automatic diagnosis strategy based on bony landmarks. • Automatic placement of the fixation plate. • Two-stage weighted multi-objective optimization based on a genetic algorithm • Novel bone protrusion evaluation considering bone contact and surfaces gaps. • Patient-specific screw optimization based on bone density information. • Capability of considering all types of common osteotomies: single-cut, opening wedge, closing wedge. Abbreviations 2D: Two-dimensional 3D: Three-Dimensional CA: Computer-assisted CAD: Computer-aided design CASPA: Planning software developed at the Computer Assisted Research and Development Group CT: Computed tomography DoF: Degrees of freedom GPU: Graphics Processing Unit GS: Gold standard OA: Optimization algorithm SSM: Statistical shape models ROM: Range of motion 1. Introduction Post-traumatic healing in non-anatomical positions (malunions) or congenital deformations of bones can cause limitations in the range of motion (ROM) of the patient, generate pain and, if not treated properly, result in severe degenerative pathologies such as osteoarthritis (Nagy et al., 2008). The current gold standard for surgical treatment of these pathologies is the restoration of the normal anatomy in a surgical procedure known as a corrective osteotomy. In this procedure, the pathological bone is cut into two or more fragments, the fragments are realigned (clinical term: reduced) to their physiological position and stabilized with an osteosynthesis implant (Schweizer et al., 2010). However, the correction of bone malunions is highly patient-specific, requiring performing a complex 6-degree-of-freedom (DoF) correction (rotation and translation) for each bone fragment, in order to restore the physiological anatomy of the patient. Therefore, corrective osteotomies are technically challenging to perform without careful diagnosis and detailed preoperative planning of the procedure. Moreover, the successful outcome of a corrective osteotomy depends also on precise intraoperative navigation of the bone reduction (Fürnstahl et al., 2016). Conventional two-dimensional (2D) preoperative planning approaches based on X-ray and Computed Tomography (CT) images fail to correctly assess the inherently three-dimensional (3D) nature of bone malunions (Schweizer et al., 2010). Another drawback of 2D preoperative planning is that it cannot be used for surgical navigation, and surgeons must rely mainly on outdated, rudimentary surgical techniques to achieve the desired bone correction. Computer-assisted 3D preoperative planning addresses those problems. Several works have established 3D preoperative planning as the state-of-the-art technique for corrective osteotomies (Athwal et al., 2003; Dobbe et al., 2011; Fürnstahl, 2010; Fürnstahl et al., 2016), due to its clear benefits in patient treatment. It allows the quantification of malunions and its corrections in all 6 DoF. 3D preoperative planning offers also the possibility of precise pre-calculation of the osteotomy fixation using 3D representations of the fixation plates and fixation screws (Dobbe et al., 2011; Miyake et al., 2012a; Schweizer et al., 2010; Schweizer et al., 2013). Lastly, it enables the translation of the preoperative planning intraoperatively by means of surgical navigation, based either on patient-specific instruments (PSI) (Fürnstahl et al., 2016; Miyake et al., 2011; Murase et al., 2008; Omori et al., 2014) or on optical navigation systems (Andress et al., 2018). Thus, the introduction of 3D preoperative planning techniques marked a paradigm shift in patient treatment, allowing performing complex surgical procedures that would not be possible using conventional techniques (Athwal et al., 2003; Dobbe et al., 2014; Fürnstahl et al., 2016; Kunz et al., 2013; Roner et al., 2017; Schweizer et al., 2013; Schweizer et al., 2016; Zdravkovic and Bilic, 1990). In our institution, the current state-of-the-art preoperative planning of long-bone osteotomies encompass the following steps: in a first step, patient-specific 3D triangular surface models (hereinafter: 3D models) of the bones are generated based on the segmentation of the CT data of the patient (Fürnstahl et al., 2008), which is part of our standard clinical procedure (Fürnstahl et al., 2008). After obtaining the patient-specific 3D bone models, a diagnosis of the malunion is performed by comparison of the pathological bone model to a reconstruction target (Figure 1A). Usually, a mirrored model of the contralateral healthy side is used as the reconstruction target. The deformity analysis allows the definition of the needed osteotomy cuts along the pathological bone model, which are simulated as shown in Figure 1B. The resulting bone fragments can then be realigned into their anatomical position by realigning the bone fragments to fit the reconstruction target (Fürnstahl, 2010; Fürnstahl et al., 2016; Murase et al., 2008). The final step of the preoperative planning procedure is the simulation of the fixation of the osteotomy by integrating 3D models of the fixation plate and screws (Figure 1C). In total, the planning process for corrective osteotomies requires the definition of 18 DoF, without including the orientation and position of the screws. This preoperative plan can be translated into the operation room by means of patient-specific navigation instruments, designed according to the preoperative plan and later 3D-printed, allowing the surgeons to perform intraoperatively step-by-step the previously simulated procedure (Fürnstahl et al., 2016; Omori et al., 2014; Rosseels et al., 2018). Despite clear advantages of 3D preoperative planning, the generation of preoperative planning solutions requires the manual calculation of the aforementioned steps by trial and error, even when using dedicated 3D planning software (Fürnstahl et al., 2016; Roner et al., 2017; Schweizer et al., 2010). The development of a clinically feasible solution also requires close collaboration between surgeons, providing the clinical knowledge, and engineers, who have the technical expertise. As the availability of surgeons needed for consultancy is often very limited, the generation of an optimal preoperative solution for extra-articular long-bone osteotomies can add up to 4 hours (Fürnstahl et al., 2016) and might involve several iterations of manual adjustments. This incurs unwanted clinical costs as a consequence of the human workload spent on manual processes. Moreover, only a reduced number of clinically feasible solutions can be investigated due to the constrained clinical timeframe. Figure 1: Overview of state-of-the-art preoperative planning of long-bone osteotomies. (A) 3D bone models obtained through segmentation of CT data. The pathological bone model is shown in white and the reconstruction target in green. Diagnosis of the malunion is done by comparison of the two bones. (B) An osteotomy plane (cyan) is created and the pathological bone model is cut, generating the proximal (white) and distal (magenta) bone fragments. (C) Bone reduction is simulated by aligning the generated fragments to the reconstruction target, and the osteotomy is fixated with a fixation plate (in blue) and its corresponding screws (in gray). One possibility to overcome those challenges is the implementation of a computer-based planning approach, able to systematically and automatically generate optimal preoperative planning solutions. The automatic generation can significantly decrease clinical costs and reduce unnecessary interaction times of the collaborators. It could also improve the quality of patient treatment by considering a larger range of parameters and solutions in order to generate better preoperative planning solutions than those obtained by a human planner. Moreover, with the current trend of the health industry towards digitalization of patient data and treatment solutions, automatic methods become an essential asset for optimal clinical treatments. However, the implementation of an automatic optimization approach is a challenging task. Early on, Zdravkovic and Bilic (1990) proposed a computer-assisted preoperative planning framework for orthopedic surgeries in order to handle the complex tasks involved in the preoperative planning process. The latter, includes multiple nonlinear, discontinuous and non-differentiable planning objectives, making their mathematical manipulation difficult. Some objectives are contradicting and tightly associated with each other, which can cause the worsening of one objective while trying to improve another one. In the field of corrective osteotomies of long bones, only a few works have tried to tackle the automatic optimization problem (Belei et al., 2007; Carrillo et al., 2017; Schkommodau et al., 2005). Although these approaches have been promising and are pioneers in the filed of automatic optimization for surgical planning of orthopedic surgeries, they still lack of solutions that can be readily applied in clinical practice without further modifications. In this paper, we present a multi-staged, multi-objective optimization approach based on the artificial intelligent methods for the generation of ready-to-use preoperative planning solutions. The system is capable of calculating solutions considering all common bone malunions, osteotomy types, and surgical approaches. We have introduce the following key contributions with respect to our previous work (Carrillo et al., 2017):  An automatic diagnosis of bone deformation based on bony landmarks, which allows automatic narrowing of the optimization search space. In our previous work, the diagnosis of the bone deformation was assumed as given.  Automatic placement of the fixation plate using information provided by statistical shape models (SSM), in contrast to the manual definition of the feasible fixation areas used in (Carrillo et al., 2017).  A novel bone protrusion evaluation that supports the generation of better solutions by considering bone contact and surfaces gaps. The use of bone protrusion represents a more realistic clinical metric than the minimization of the cut surface used in our previous work.  In contrast to the heuristic approach employed in our previous work, in this paper, we have developed a patient-specific screw purchase optimization based on bone density information extracted from CT data.  The capability of considering all types of common osteotomies (single-cut, opening wedge, closing wedge) along the entire Radius, in order to generate solutions that would be difficult to achieve for a human planner within a reasonable time. The algorithm presented in (Carrillo et al., 2017) was only able to handle distal radius osteotomies and was not capable of generating single-cut solutions.  We have also reduced the computational effort of the strategy presented in (Carrillo et al., 2017) by applying a different multi-stage approach.  Finally, we have provided a profound clinical and mathematical explanation of all the optimization objectives, which was previously missing. The developed optimization framework was validated clinically against state-of-the-art preoperative planning solutions on a consecutive series of patients with forearm malunions, previously treated at our institution. In the following, we will give a brief overview of existing preoperative planning approaches for orthopedic surgeries. In section 2, we describe in detailed the proposed optimization framework. In section 3, dataset, experimental set-up, and results are presented. In section 4 a discussion about results, impact, and limitations of this work is given. Finally, in section 5, we draw final conclusions and give an outlook about future work. 1.1. Related Work The first step towards the generation of a complete 3D preoperative planning is the patient diagnosis done through the analysis of the bone malunion. In 3D preoperative planning, patient-specific 3D bone models generated from the segmentation of multi-planar data (Lorensen and Cline, 1987) are superimposed over a healthy reconstruction template using semi-automatic registration methods (Kawakami et al., 2002; Schenk et al., 2016; Schweizer et al., 2010; Vlachopoulos et al., 2017). Afterward, the bone deformation is quantified by means of clinical and mathematical measurements (Fürnstahl et al., 2016; Gosse et al., 1997; Murase et al., 2008; Subburaj et al., 2010), providing the basis for the preoperative planning of the surgical correction. Here, 4 main objectives have to be considered: osteotomy cut plane, reduction of the bone fragments, position of the fixation plate and direction and position of fixation screws. Current state-of-the-art of 3D preoperative planning for corrective osteotomies uses dedicated planning software for the manual calculation of each of these objectives. We refer the interested readers to (Fürnstahl, 2010; Fürnstahl et al., 2016; Schweizer et al., 2016; Vlachopoulos et al., 2015) for more information about the 3D planning tool. The latter facilitates the process of generation of a 3D preoperative plan, however the basic primitives operations needed for creating a preoperative plan are cutting Boolean operations and the possibility for interactive transformation, which would be also available in any dedicated CAD software. The process of creating each of the steps involved in the preoperative plan is difficult to control even by skilled engineers, as any change done to one of the parameters (e.g., position of fixation screws, inclination of osteotomy plane) must be manually propagated across the different objectives (Athwal et al., 2003; Belei et al., 2007; Bilic et al., 1994; Carrillo et al., 2017; Fürnstahl, 2010). Current-state-of-the-art planning (Fürnstahl et al., 2016; Roner et al., 2017; Vlachopoulos et al., 2015) is also incapable of handling contradicting objectives, meaning for example that an improvement in the accuracy of the bone reduction can subsequently deteriorate the position of the fixation plate or generate solutions with non-feasible osteotomy cuts. Similarly, an improvement in the position of the osteotomy cut can cause a less fitting position for the placement of the fixation plate. Existing automatic methods have only solved a reduced problem set, failing to include all objectives needed for a ready-to-use clinical solution , i.e., osteotomy cut plane, reduction of the bone fragments, position of the fixation plate and direction and position of fixation screws. Such is the case of previously described methods (Eck et al., 1990; Menetrey and Paul, 2004), where an automatic calculation was performed only for the osteotomy plane. Eck et al. (1990) calculated the osteotomy plane in femoral head reduction planning using a nonlinear optimization algorithm, based on a least-square-approximation solver. Menetrey and Paul (2004) did a similar parametrization of the osteotomy plane and wedge size for osteotomies around the knee, optimizing only the reduction alignment. Schkommodau et al. (2005) developed a multi-objective optimization strategy for corrective osteotomies of lower extremities that considered the following optimization objectives: leg length, translation, antetorsion, and angulation. Their method solved a simplified osteotomy sub-problem with only 12 DoF, which did not include the position of the fixation plate or the fixation screws into the optimization process. The multi-objective problem was solved by a sequential quadratic programming algorithm and considered the influence of all these objectives in the preoperative planning. The approach was later extended by Belei et al. (2007) to account for different osteotomy types (closing wedge, opening wedge, single cut) and to consider also double osteotomy solutions. These approaches represent the first attempt to automate the planning of corrective osteotomy. However, the planning was based on simplified geometry rather than on patient anatomy and it relied on intraoperatively calibrated fluoroscopic datasets. Learning-based methods have proven to be effective in similar applications of medical image processing techniques (Criminisi and Shotton, 2013; Esfandiari et al., 2018; Tschannen et al., 2016). In the field of shoulder arthroplasty, Tschannen et al. (2016) presented an automatic algorithm for preoperative planning of the resection plane for arthroplasty of the proximal humeral head, based on random regression forests. The approach allowed controlling the orientation, position, and size of the prosthetic humeral head in relation to the humeral shaft, using a direct mapping between the CT image and the parameters of the resection plane. The estimation of the plane position using CT data could be of interest to speed-up the convergence of optimization algorithms. Another interesting application is found in the field of spine surgery, where the problem of pedicle screw placement and pose estimation has been extensively studied (Farshad et al., 2017; Scheufler et al., 2011). Also, Esfandiari et al. (2018) proposed an algorithm based on convolutional neural networks for the 6-DoF estimation of the screw position and direction using intraoperative fluoroscopy data and estimation of the bone density. Similar approaches, taking into account the bone density information of the patient, should be considered in long-bone osteotomies to increase the quality of the osteotomy fixation. To the best of our knowledge, only our previous work (Carrillo et al., 2017) deals with an 18-DoF optimization problem for the generation of preoperative planning solutions of corrective osteotomies. The optimization approach presented in (Carrillo et al., 2017) has been taken as the basis for the core optimization algorithm (section 2.3.3) of the herein presented framework. 2. Methods We have developed an optimization framework for the generation of ready-to-use preoperative planning solutions for corrective osteotomies. A complete overview of our approach is given in Figure 2. The framework receives the 3D bone models of the patient, the reconstruction target and the fixation plate as an input (section 2.1). Afterward, an automatic diagnosis of the malunion is performed by identification of the pathological area, encoding also feasible plate regions and associated clinical constraints (section 2.2). This information is used by the optimization module, where the generation and optimization of preoperative planning solutions take place (i.e., the simulation of the osteotomy cut, realignment of the fragment and virtual fixation of the osteotomy) using a multi-stage optimization strategy (section 2.3). A genetic algorithm approach is used for optimization in which the objectives are encoded in a real-valued chromosome. The optimization is controlled by means of fitness evaluation functions, which provide quality measurements for the attainment of the different objectives. The output of the framework is a complete preoperative planning solution, including position and orientation of the osteotomy plane, the 6-DoF transformation to achieve the reduction, and the position and orientation of the fixation plate and screws. Our approach has been designed to cover all long-bone malunions, however, we have focused the implementation of our method on the radius due to the high rate of radius osteotomy procedures performed yearly at our institution. Figure 2: Overview of the automatic optimization framework for the 3D planning of long bone corrective osteotomy surgeries. 2.1. Input data generation In our institution, patient-specific 3D models of bone anatomy are used in the standard treatment workflow for complex malunions. For a better understanding, we briefly describe the model generation process although it is not an integral part of our optimization framework. Patient-specific 3D bone models of the forearm were generated from CT data (slice thickness 1 mm; 120 kV; Philips Brilliance 64 CT, Philips Healthcare, The Netherlands) using thresholding and region-growing algorithms of commercial segmentation software (Mimics Medical, Version 19.0, Materialise 2016, Leuven, Belgium). 3D models are generated using the Marching Cube algorithm (Lorensen and Cline, 1987) and given in the form of triangular surface meshes (stereolithographic models; STL) as described by (Roscoe, 1988). In our institution, the mirrored model of the contralateral bone is always used as the anatomical (healthy) reconstruction target. As the last input parameter, a feasible osteosynthesis implant (fixation plate) must be chosen, used for the fixation of the bone fragments after osteotomy. All 3D models were transformed into a common anatomical reference frame. In the case of the radius, we have used the anatomical coordinate system described in (Vlachopoulos et al., 2015). The coordinate system is � � � � ⃑ 𝐶𝐶𝐶𝐶 and it is oriented as shown in Figure 3A. Different anatomical and clinically relevant landmarks denoted by regions were defined and used throughout the entire pipeline. 7 point sets of landmarks regions, of 5 mm radius size, were automatically generated on the distal parts of both, pathological radius {𝐿𝐿 𝑃𝑃 , … ,𝐿𝐿 𝑃𝑃 } and reconstruction 1 7 target {𝐿𝐿 𝑅𝑅 , … ,𝐿𝐿 𝑅𝑅 } as depicted on Figure 3B. In our implementation, each set 𝐿𝐿 𝑃𝑃 and 𝐿𝐿 𝑅𝑅 contains a total of 1 7 𝑙𝑙 𝑙𝑙 𝑁𝑁 = 50 points, and each point set 𝐿𝐿 𝑃𝑃 is in correspondence to the reciprocal point set 𝐿𝐿 𝑅𝑅 . 𝑙𝑙 𝑟𝑟 � � � � � ⃗ Figure 3: Anatomical Axis and Bony Landmarks. (A) Anatomical coordinate system; coordinate axis 𝑪𝑪𝑪𝑪 is depicted by 𝒙𝒙 � � � � � ⃗ the red arrow and represents the radioulnar direction, the coordinate axis 𝑪𝑪𝑪𝑪 is depicted by the green arrow and represents 𝒚𝒚 � � � � � ⃗ the axial direction (long bone axis), and coordinate axis 𝑪𝑪𝑪𝑪 is shown in blue and indicates the dorsal-palmar direction. (B) 𝒛𝒛 Anatomical landmarks of the distal radius: (i) dorsal distal edge of the sigmoid notch; (ii) palmar distal edge of the sigmoid notch; (iii) center of distal radius styloid; (iv) center of scaphoid facet; (v) center of lunate facet; (vi) Lister’s tubercle and (vii) center of radius’ watershed line. 2.2. Automatic Diagnosis of Bone Malunion A common clinical practice for the 3D diagnosis of a bone malunion is to align the healthy proximal and distal joints of the pathological bone to the corresponding regions of the reconstruction target (Schweizer et al., 2010b). This process allows obtaining a visual approximation of the deformed area of the pathological bone and it defines also the most proximal and most distal regions of the bone at which an osteotomy cut can be performed. The bone malunion can then be quantified in 3D by the relative 4𝑥𝑥 4 transformation matrix between the proximally and distally aligned bone models (Figure 4). The transformation matrix encodes the translation � � � � ⃑ and rotational deviation of the malunion in all 6 DoF with respect to the anatomical coordinate system 𝐶𝐶𝐶𝐶 (Schweizer et al., 2010; Vlachopoulos et al., 2017). Based on this well-established principle, we have developed an automatic method for the calculation of the deformity region, the deformity profile, and the quantification of the bone malunion. Figure 4: Method for automatic diagnosis of malunions. (A-B) The pathological bone P is proximally (A) and distally (B) 𝒑𝒑 aligned to the reconstruction target R. (C) Example calculation of windows function for 𝑷𝑷 ; function 𝑫𝑫𝒑𝒑 is calculated along � � � � � ⃗ ( ) the 𝑪𝑪𝑪𝑪 of the bone, from 𝑳𝑳 to 𝑳𝑳 , with window function 𝑾𝑾 (∙) of thickness 𝒕𝒕 = 𝟏𝟏 𝒎𝒎𝒎𝒎 ,∞,∞ . (D) Deformity 𝒚𝒚 𝒙𝒙 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 curves; Lower envelopes 𝒍𝒍𝒍𝒍𝒍𝒍 (bold black signal) of 𝑫𝑫𝒑𝒑 (𝑷𝑷 ) (blue signal), and 𝒍𝒍𝒍𝒍𝒍𝒍 (bold magenta signal) of 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 𝑫𝑫𝒑𝒑 (𝑷𝑷 ) (gray signal), are calculated, corresponding to both proximal and distal alignments. (E) Deformity region is shown in orange, between the most proximal and most distal planes obtained from the deformity threshold 𝒅𝒅 applied to 𝑻𝑻𝑻𝑻 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 𝒅𝒅 𝒎𝒎 𝒕𝒕𝒅𝒅 𝒍𝒍𝒍𝒍𝒍𝒍 and 𝒍𝒍𝒍𝒍𝒍𝒍 . In the automatic diagnosis, the pathological bone model 𝑷𝑷 is coarsely aligned to the reconstruction target 𝑹𝑹 using a landmark-based registration method, where each 𝐿𝐿 𝑃𝑃 of the pathological bone is registered to its 𝑙𝑙 𝒎𝒎𝒎𝒎 𝒎𝒎𝒎𝒎𝒎𝒎 reciprocal 𝐿𝐿 𝑅𝑅 of the reconstruction target, using standard iterative closest point (ICP) (Besl and McKay, 1992). 𝑙𝑙 Afterwards, the minimum and maximum points 𝐿𝐿 and 𝐿𝐿 are respectively calculated, on the bone length 𝑚𝑚 𝑚𝑚 � � � � ⃗ axis 𝐶𝐶𝐶𝐶 with respect to 𝑃𝑃 . Let 𝑦𝑦 𝛿𝛿 (𝑥𝑥 ,𝑌𝑌 ) = arg min� 𝑥𝑥 − 𝑦𝑦 � (Eq. 1) 𝑗𝑗 𝑦𝑦 ∈ 𝑌𝑌 𝑗𝑗 be the nearest neighbor function yielding the point 𝑦𝑦 of a point set or line 𝑌𝑌 with the smallest Euclidean distance 𝑗𝑗 ( ) to a point 𝑥𝑥 . 𝐿𝐿 and 𝐿𝐿 , and the line segment L 𝜆𝜆 between them can then be defined as 𝑚𝑚 𝑚𝑚 � � � � ⃗ � � � � ⃗ � � � � ⃗ � � � � ⃗ 𝐿𝐿 = 𝛿𝛿 � ,𝐶𝐶𝐶𝐶𝑝𝑝 � | 𝑝𝑝 = min� ∙ 𝐶𝐶𝐶𝐶 𝑝𝑝 � , 𝐿𝐿 = 𝛿𝛿 � ,𝐶𝐶𝐶𝐶𝑝𝑝 � | 𝑝𝑝 = max� ∙ 𝐶𝐶𝐶𝐶 𝑝𝑝 � , 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑝𝑝 ∈ 𝑃𝑃 𝑝𝑝 ∈ 𝑃𝑃 𝑖𝑖 𝑖𝑖 and 𝐿𝐿 (𝜆𝜆 ) = 𝐿𝐿 + 𝜆𝜆 (𝐿𝐿 −𝐿𝐿 ) | 0 ≤ 𝜆𝜆 ≤ 1, 𝑚𝑚 𝑚𝑚 𝑚𝑚 where · denotes the scalar product. Following the recommendation of (Vlachopoulos et al., 2017) for the size of the proximal and distal alignment regions, L(𝜆𝜆 ) is used to determine the points of 𝑃𝑃 that lie in the most proximal and distal 30% with respect to � � � � ⃗ 𝐶𝐶𝐶𝐶 : 𝑦𝑦 � � � � ⃗ � � � � ⃗ � � � � ⃗ { 𝑝𝑝 | 𝑝𝑝 ∙ 𝐶𝐶𝐶𝐶 ≤ 𝐿𝐿 (0.3)∙ 𝐶𝐶𝐶𝐶 } and { 𝑝𝑝 | ||𝑝𝑝 ∙ 𝐶𝐶𝐶𝐶 − 𝐿𝐿 || ≤ ||𝐿𝐿 (0.3) − 𝐿𝐿 ||} 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑦𝑦 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑚𝑚 𝑚𝑚 These points are then used to align the proximal and distal joints of 𝑷𝑷 to 𝑹𝑹 (Figure 4A). The alignment is done using ICP registration (Besl and McKay, 1992) and the proximally and distally aligned models are denoted by 𝒑𝒑 𝒅𝒅 𝑷𝑷 and 𝑷𝑷 , respectively. 𝑝𝑝 𝑑𝑑𝑝𝑝 Afterward, the deformity profiles 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) (Figure 4D, blue signal) and 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) (Figure 4D, grey signal) are obtained from the deformity function 𝐷𝐷 𝑝𝑝 (𝑃𝑃𝐶𝐶 ) = �𝑊𝑊 (𝑙𝑙 (𝜆𝜆 ),𝑝𝑝 ) − 𝛿𝛿 � (𝑙𝑙 (𝜆𝜆 𝑊𝑊 ),𝑝𝑝 ),𝑊𝑊 (𝑙𝑙 (𝜆𝜆 ),𝑅𝑅 )� , (Eq. 2) 𝑚𝑚 𝑚𝑚 � | | 𝑃𝑃𝐶𝐶 | | 𝑚𝑚 =1,.., 𝑃𝑃𝑃𝑃 𝜆𝜆 =0,0.1,..,1 | | ( ) where 𝑃𝑃𝐶𝐶 is a pointset, 𝑃𝑃𝐶𝐶 is the cardinality of the pointset and 𝑊𝑊 · denotes the window function { | 𝑝𝑝 ∈ 𝑃𝑃𝐶𝐶 𝑐𝑐 −𝑡𝑡 < 𝑝𝑝 < 𝑐𝑐 +𝑡𝑡 } 𝑚𝑚 𝑚𝑚 ( ) 𝑊𝑊 𝑐𝑐 ,𝑃𝑃𝐶𝐶 = � . (Eq. 3) ∅ ∶ 𝑜𝑜𝑡𝑡ℎ 𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 𝑐𝑐 is the center of the window and 𝑡𝑡 is an adjustable user-defined thickness. As shown in Figure 4C, we have � � � � ⃗ used a thickness 𝑡𝑡 = 0.01 𝑚𝑚𝑚𝑚 , and we have taken 𝑡𝑡 and 𝑡𝑡 to be large enough to enclose 𝑷𝑷 and 𝑹𝑹 along 𝐶𝐶𝐶𝐶 𝑦𝑦 𝑚𝑚 𝑧𝑧 𝑚𝑚 � � � � ⃗ and 𝐶𝐶𝐶𝐶 . 𝑧𝑧 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 In order to reduce noise in the signals, we have implemented filtering using the lower envelopes 𝑜𝑜𝑒𝑒𝑙𝑙 (Figure 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑝𝑝 𝑑𝑑 4D, bold black line) and 𝑜𝑜𝑒𝑒𝑙𝑙 (Figure 4D, bold magenta line) of 𝐷𝐷 𝑝𝑝 (𝑃𝑃 ) and (𝑃𝑃 ) , respectively. Based on 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑜𝑜𝑒𝑒𝑙𝑙 and 𝑜𝑜𝑒𝑒𝑙𝑙 , a deformity threshold 𝑑𝑑 = 10% defines the part of the bone in which the deformity profile 𝑇𝑇 ℎ deviates more than 10% from their minimum value (Figure 4E). This region defines the most proximal and most distal sites for feasible osteotomy cuts. Finally, the malunion is quantified by a 4𝑥𝑥 4 homogenous transformation matrix 𝑀𝑀 representing the 𝑙𝑙 difference between the proximally and distally aligned 𝑷𝑷 . 𝑀𝑀 contains the information about the rotational and 𝑙𝑙 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 � � � � ⃗ � � � � ⃗ translational error of the bone malunion with respect to 𝐶𝐶𝐶𝐶 . The rotational deviation along each direction of 𝐶𝐶𝐶𝐶 is obtained by applying the Euler-Rodrigues formula (Schneider and Eberly, 2002) to the rotational component of 𝑀𝑀 , which yields the Euler angles 𝜑𝜑 ,𝜑𝜑 and 𝜑𝜑 . The translational deviation along each coordinate axis is 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑧𝑧 obtained from the translational parts 𝑇𝑇 ,𝑇𝑇 and 𝑇𝑇 of 𝑀𝑀 (Schneider and Eberly, 2002). Thus, the 6-DoF 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑙𝑙 � ⃗ quantification of bone malunion is given by the vector 𝑉𝑉 = � 𝜑𝜑 ,𝜑𝜑 ,𝜑𝜑 ,𝑇𝑇 ,𝑇𝑇 ,𝑇𝑇 � . 𝑚𝑚𝑚𝑚 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑚𝑚 𝑦𝑦 𝑧𝑧 2.3. Multi-stage optimization Our optimization approach comprises three parts. First, each clinical goal is characterized mathematically by means of optimization parameters, boundaries, and constraints. These parameters are used for the automatic generation of solutions at each optimization step, and the boundaries and constraints are used to control the search space of the algorithm. Our mathematical translation of the clinical goals into optimization objectives is explained in section 2.3.1. Secondly, the optimization algorithm requires a way to evaluate the performance of the optimization objectives. In our algorithm, this evaluation is performed by tailored fitness functions, developed for the minimization or maximization of the optimization parameters. The appropriate definition of the fitness function is crucial for the generation of feasible solutions and convergence of the algorithm. We give a detailed description of the fitness functions in section 2.3.2. Lastly, an optimization strategy must be chosen for the generation of solutions. In our case, the optimization algorithm had to be able to handle multiple objectives with 18 DoF and nonlinear constraints. We have implemented a multi-stage optimization approach based on the multi-objective genetic algorithm NSGA-II (Deb et al., 2000), further explained in section 2.3.3. 2.3.1. Optimization Objectives A complete preoperative solution for long bone osteotomies entails the calculation of 4 clinical objectives, as shown in Figure 5: (A) The position and orientation of the osteotomy cut plane (6 DoF), (B) the translation and rotation of the generated bone fragments for the reduction to the reconstruction target (6 DoF), (C) the position and orientation of the fixation plate (6 DoF) and (D) the position and orientation of the plate’s screws into the reduced fragments. The definition of each of these objectives requires more than a simple mathematical translation into adjustable parameters. We have to carefully decide on the parameters to control such that the parameter space can be minimal, and the controllability of the objectives can be directly evidence from a parameter change. In this section, we describe the strategy on the mathematical parameters, constraints, and boundaries associated with each objective. 𝑚𝑚𝑚𝑚 𝑚𝑚𝑚𝑚 Figure 5: Optimization objectives. (A) Position and orientation of the osteotomy cut plane (in cyan) along the pathological bone (in white). (B) Reduction alignment of the generated bone fragments (in white) with respect to the reconstruction target (in green). (C) Position of the fixation plate (in gray) with respect to the proximal (in white) and transformed distal fragment (in magenta). (D) Position and orientation of the screws of the fixation plate (in blue and gray), used for the stabilization and fixation of the bone fragments. A. Osteotomy Plane (Figure 5A). The osteotomy plane determines the size and location of the osteotomy and the orientation of the bone cut. It is represented by the position 𝑙𝑙𝑝𝑝 (3 DoF) and normal 𝑛𝑛� ⃗ (3 DoF), with a 𝑝𝑝𝑙𝑙 � � � � ⃗ total of 6 DoF. Regarding 𝑙𝑙𝑝𝑝 , we optimize only the position 𝑙𝑙𝑝𝑝 relative to the long bone axis 𝐶𝐶𝐶𝐶 (Table 1). 𝑦𝑦 𝑦𝑦 � � � � ⃗ In extra-articular osteotomies, only the translation of the osteotomy plane along the bone length axis 𝐶𝐶𝐶𝐶 has 𝑦𝑦 � � � � ⃗ � � � � ⃗ an effect on the osteotomy, thus the translation of the osteotomy plane with respect to 𝐶𝐶𝐶𝐶 and 𝐶𝐶𝐶𝐶 is 𝑚𝑚 𝑧𝑧 redundant. Together with the reduction alignment, the osteotomy plane determines the type of osteotomy and the shape of the wedge between the generated fragments. In our approach, we have considered all common osteotomy types (Figure 6): opening-wedge, closing wedge and single-cut osteotomy. Figure 6A shows the situation of the bone previous to realignment and after performing the osteotomy cut (in cyan). In an opening-wedge osteotomy (Figure 6B) a bone cut is made and the generated bone fragments are reduced to their anatomical position, yielding a wedge-shaped opening. A too large gap may prevent healing or may result in implant failure due to instability. To avoid this, large gaps are filled with structural bone graft in the surgery to support healing (Fernandez, 1982; Fürnstahl et al., 2016; Murase et al., 2008). Oppositely, a closing-wedge osteotomy (Figure 6C) refers to an osteotomy where a bone wedge has to be removed to complete the bone reduction. This may occur due to an overlap between the two generated bone fragments after their reduction to the reconstruction target (Fürnstahl et al., 2016; Murase et al., 2008). Our algorithm calculates first the necessary reduction of the bone fragments and afterward the overlap is defined by simulating a second planar cut. Additionally, our approach has also the feature of calculating single-cut osteotomies (Figure 6D) (Dobbe et al., 2017; Sangeorzan et al., 1989). A single cut osteotomy is always preferred against closing wedge solutions due to better bone contact and faster bone healing (Merle d’Aubigné and Descamps, 1952; Roner et al., 2017; Sangeorzan et al., 1989). This type of osteotomy is technically more complex to plan manually, as the bone reduction is achieved by sliding and rotating the generated bone fragments along the osteotomy cut plane (Sangeorzan et al., 1989). We achieve the support of single-cut osteotomies in our framework by constraining the closing wedge osteotomy planes to be parallel � � � � ⃗ and overlying. In this case, the wedge is limited to a maximum overlap of 1 mm along the 𝐶𝐶𝐶𝐶 . 𝑦𝑦 In each optimization stage, we use the polygon clipping algorithm of (Vatti, 1992) to cut model P by plane 𝒑𝒑𝒑𝒑 𝒅𝒅 𝒅𝒅𝒎𝒎 𝒕𝒕 (𝑙𝑙𝑝𝑝 ,𝑛𝑛 ), yielding potential candidates of the proximal and distal fragments 𝑷𝑷 and 𝑷𝑷 , respectively. 𝑝𝑝𝑙𝑙 The type of osteotomy is controlled by the maximum and minimum allowed wedge size using nonlinear constraints. Figure 6: Osteotomy Types. (A) Situation after performing the osteotomy cut (shown in cyan) and before reduction to the reconstruction target. The proximal fragment is shown in white and the distal fragment in magenta. (B) Opening wedge: distal fragment (magenta) has been reduced to the reconstruction target, generating an opening between both fragments. (C) Closing wedge: in order to complete the osteotomy, after simulated reduction of the distal fragment (magenta), a second planar cut must be done to remove the wedge of the overlapping bone region (shown in red). (D) Single cut: reduction of the distal fragment is achieved by fragment rotation and translation along the osteotomy plane. 𝒅𝒅 𝒅𝒅𝒎𝒎 𝒕𝒕 B. Reduction Alignment (Figure 5B). Once the osteotomy plane cut is executed, 𝑷𝑷 must be rotated and translated to match the given reconstruction target. This process of bone reduction represents one of the most important objectives for clinicians, because the precision of the reduction significantly influences the surgical outcome. Mathematically, the reduction alignment can be represented by a 4x4 homogenous transformation � ⃗ matrix (A ), constructed from the rotation 𝑅𝑅 = � ,𝜑𝜑 𝜑𝜑 ,𝜑𝜑 � and translation 𝑇𝑇 = � ,𝑇𝑇 𝑇𝑇 ,𝑇𝑇 � parameters f 𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑓𝑓 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 along each axis (6 DoF; Table 1). C. Position of Fixation Plate (Figure 5C). Another factor influencing the success of the treatment is the position of the fixation plate. A poorly fitting plate results in a less controlled procedure, decreases stability of the reduction and even leads to inaccurate reduction, delayed bone healing or failure of the fixation plate (Blecha et al., 2005; Fürnstahl et al., 2016). Mathematically, the position of the fixation plate along the bone entails a 6 DoF task, which can be encoded in the optimization by a 4x4 homogenous transformation matrix 𝒍𝒍𝒙𝒙 � ⃗ 𝐴𝐴 . The latter is constructed from the rotation 𝑅𝑅 = � ,𝜃𝜃 ,𝜃𝜃 𝜃𝜃 � and t he translation 𝑇𝑇 = � ,𝑇𝑇 𝑇𝑇 ,𝑇𝑇 � of 𝑝𝑝 𝑝𝑝 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑝𝑝 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 � � � � ⃗ the fixation plate (Table 1), relative to the reference coordinate system 𝐶𝐶𝐶𝐶 . D. Screw Position and Orientation (Figure 5D). The quality of the solution is also considerably dependent on the position and orientation of the screws, which are required to achieve a stable bone-plate interface. Several studies have shown how the length of the fixation screw can affect the stability of the osteotomy and how poorly placed screws could lead to soft-tissue inflammation and ligament tearing (Blecha et al., 2005; Dougherty et al., 2008; Reitman et al., 2004). Mathematically, the position and orientation of a fixation screw S is directly related to the position and orientation of the fixation plate, and can be described by � ⃗ S = h A � + R t � (Eq. 4) i s p s s i i i � ⃗ where h is the screw axis, A is the 4𝑥𝑥 4 homogenous transformation matrix of the fixation plate position, s p R is a 4𝑥𝑥 4 rotation matrix encoding orientation of the screw relative to the plate and t is a constant s s i i translation vector containing the relative position of the screw hole with respect to the plate (Table 1). For uniaxial locking screw plate systems, where the orientation of the screw is already given by the manufacturer, � � � � � ⃗ � � � � � ⃗ R is assumed to be constant, i.e., R = dır I , where dır is a constant direction vector. In the s s s 4x4 s i i i i � ⃗ optimization, the orientation of each screw is described by a vector θ , whose components are the Euler � ⃗ angles derived from R . The number of vectors θ is determined by the number of screws 𝑁𝑁 = |{𝑒𝑒 }| of the s s 𝑑𝑑𝑚𝑚 𝑚𝑚 i i fixation plate. Once the optimization objectives have been defined, the optimization parameters are encoded directly as decimal values into a parametrization vector. The vector is further denoted as chromosome 𝑥𝑥⃗, defined as � ⃗ ⃗ ⃗ 𝑥𝑥⃗ = � 𝑛𝑛 𝑙𝑙𝑝𝑝 𝑛𝑛 𝑛𝑛 𝜑𝜑 𝜑𝜑 𝜑𝜑 𝑇𝑇 𝑇𝑇 𝑇𝑇 𝜃𝜃 𝜃𝜃 𝜃𝜃 𝑇𝑇 𝑇𝑇 𝑇𝑇 θ 𝜃𝜃 … 𝜃𝜃 � , (Eq. 5) 𝑦𝑦 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑚𝑚 𝑦𝑦 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 s 𝑑𝑑 𝑑𝑑 𝑧𝑧 𝑧𝑧 𝑥𝑥 𝑦𝑦 𝑧𝑧 1 2 𝑁𝑁 𝑠𝑠𝑖𝑖 which will be used by the fitness functions to control changes in the optimization objectives. Table 1: Summary of optimization objectives and their associated fitness functions, optimization parameters and constraints Optimization Fitness Optimization Challenges Parameters Constraints Objective Function Parameters Represented by 𝑨𝑨 ; 𝒇𝒇 Landmark-  Accuracy of joint  Rotation 𝑅𝑅 𝑓𝑓 based Reduction surface Maximal allowed Registration (𝜑𝜑 , 𝜑𝜑 , 𝜑𝜑 ) 6 𝑚𝑚 𝑦𝑦 𝑧𝑧 Alignment  Dependent on reduction error Error � � ⃗  Translation T landmark weighting (𝑓𝑓 ) (𝑇𝑇 , 𝑇𝑇 , 𝑇𝑇 ) 𝑓𝑓𝑚𝑚 𝑓𝑓𝑦𝑦 𝑓𝑓𝑧𝑧  Avoid longitudinal  For single-cut and intraarticular 𝑟𝑟𝑃𝑃 𝑃𝑃 𝑑𝑑𝑝𝑝𝑑𝑑𝑃𝑃  Position 𝑝𝑝𝑙𝑙  𝑛𝑛� ⃗ ∥ 𝑛𝑛� ⃗ 𝑝𝑝𝑙𝑙 𝑝𝑝𝑙𝑙 cuts (𝑝𝑝𝑙𝑙 ,𝑝𝑝𝑙𝑙 , 𝑝𝑝𝑙𝑙 ) Osteotomy Bone Protrusion 𝑚𝑚 𝑦𝑦 𝑧𝑧 6  Size of bone wedge  Osteotomy type and  Normal 𝑛𝑛� ⃗ Plane (𝑓𝑓 ) 𝑝𝑝 𝑙𝑙  Range given by location malunion (𝑛𝑛 , 𝑛𝑛 , 𝑛𝑛 ) diagnosis of section pl 𝑝𝑝 𝑙𝑙 𝑝𝑝 𝑙𝑙 𝑥𝑥 𝑦𝑦 𝑧𝑧 influences reduction 2.2 alignment  Bone plate Represented by 𝑨𝑨 ; 𝒑𝒑  Gaps / steps between penetration Distance  Rotation 𝑅𝑅 𝑝𝑝 bone fragments and  Maximal allowed Fixation Plate - Position (𝜃𝜃 , 𝜃𝜃 , 𝜃𝜃 ) plate to be avoid 𝑚𝑚 𝑦𝑦 𝑧𝑧 6 distance to bone Fixation Plate Bone Fragments  Stable alignment � ⃗  Osteotomy site:  Translation 𝑇𝑇 𝑝𝑝 (𝑓𝑓 )  Various plate types palmar, dorsal, (𝑇𝑇 , 𝑇𝑇 , 𝑇𝑇 ) 𝑝𝑝 𝑚𝑚 𝑝𝑝 𝑦𝑦 𝑝𝑝 𝑧𝑧 lateral.  Proximity to joint � ⃗ ⃗ ⃗ � θ 𝜃𝜃 … 𝜃𝜃 � s 𝑑𝑑 𝑑𝑑  Minimal distance to 1 2 𝑁𝑁 𝑖𝑖𝑠𝑠 area 3 x 𝑁𝑁 𝑚𝑚𝑑𝑑 osteotomy plane Position and  Avoid penetration (𝑁𝑁 is the 𝑚𝑚𝑑𝑑 𝜃𝜃 = (𝜃𝜃 ,𝜃𝜃 ,𝜃𝜃 ) Screw Purchase 𝑑𝑑 𝑑𝑑 𝑑𝑑 𝑑𝑑  Number of screws 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑚𝑚 𝑦𝑦 𝑧𝑧 with osteotomy Orientation of number of (𝑓𝑓 ) inside the bone 4 for 𝑒𝑒 = 1, … ,𝑁𝑁 𝑚𝑚𝑑𝑑 Screws plane fixation  Bi-cortical solutions w.r.t plate coordinate  Patient-specific bone screws) preferred system and 𝐴𝐴 𝑝𝑝 density 2.3.2. Fitness functions and constraints Genetic algorithms employ minimization/maximization functions for the quantitative evaluation of the optimization parameters 𝑥𝑥⃗ during the optimization process. These functions are known as fitness functions and their definition is always challenging as they must be capable of providing a quantitative measurement of the quality of the different optimization objectives. The design of these fitness functions is a very challenging task, as we must consider not only mathematically solid functions that can be used by an optimization algorithm, but we must also guarantee that these functions reflect the clinical needs. For example, we cannot measure the quality of the osteotomy cut by simply measuring the inclination angle of the cut, because the bone contact and the bone protrusion of the generated fragments are of more clinical relevance than the orientation of the plane. In this section, we explain the design of the 4 fitness functions that we have developed based on clinically relevant measurements. Control of bone reduction alignment From the clinical point of view, some bone regions (i.e., joint regions) are more important to be precisely aligned to the reconstruction target than others. Moreover, in some cases, some bone regions inevitably deviate due to the pathology. We have developed a landmark-based error measurement to fine-control the reduction alignment by landmark regions 𝐿𝐿 𝑃𝑃 and 𝐿𝐿 𝑅𝑅 of the pathological bone and the reconstruction, respectively (Figure 3B). A 𝑙𝑙 𝑙𝑙 weight 𝑒𝑒 between 0 and 1 was assigned to all points of 𝐿𝐿 𝑃𝑃 , with ∑𝑒𝑒 = 1. Landmark regions located on the 𝑙𝑙 𝑙𝑙 𝑙𝑙 joint surface (i-v; Figure 3B) were given a bigger weight due to the importance of the joint for the reconstruction { } accuracy, which yields the following weight distribution 𝑒𝑒 , … ,𝑒𝑒 = 0.18 , 𝑒𝑒 ,𝑒𝑒 = 0.05 . In this way, the 1 5 6 7 quality of the reduction alignment is controlled by |{𝐿𝐿𝑃𝑃 }| 𝑁𝑁 =50 1 1 2 𝑓𝑓 = � 𝑒𝑒 � � � 𝑝𝑝 𝐴𝐴− 𝑞𝑞 � 𝑒𝑒 ℎ𝑒𝑒𝑒𝑒𝑒𝑒 𝑝𝑝 ∈ 𝐿𝐿 𝑃𝑃 ,𝑞𝑞 ∈ 𝐿𝐿 𝑅𝑅 (Eq. 6) 1 𝑙𝑙 𝑓𝑓 𝑚𝑚 𝑚𝑚 𝑚𝑚 𝑙𝑙 𝑚𝑚 𝑙𝑙 |{𝐿𝐿𝑃𝑃 }| 𝑁𝑁 𝑚𝑚 =1 𝑙𝑙 =1 which is the average weighted 𝑅𝑅𝐶𝐶𝑀𝑀𝑅𝑅 among the Euclidean distances between all 𝐴𝐴 -transformed points 𝑝𝑝 of the 𝑓𝑓 𝑚𝑚 pathological landmark regions and the points 𝑞𝑞 of the landmark regions of the reconstruction target. |{𝐿𝐿𝑃𝑃 }| is 𝑚𝑚 the cardinality of the set 𝐿𝐿𝑃𝑃 . The points are in correspondence, meaning that 𝑝𝑝 and 𝑞𝑞 have the same position 𝑚𝑚 𝑚𝑚 relative to their landmark centers. 𝐴𝐴 is constructed from 𝑅𝑅 = (𝜑𝜑 ,𝜑𝜑 ,𝜑𝜑 ) and 𝑡𝑡 = (𝑡𝑡 , 𝑡𝑡 ,𝑡𝑡 ) obtained 𝑓𝑓 𝑓𝑓 𝑚𝑚 𝑦𝑦 𝑧𝑧 𝑓𝑓 𝑚𝑚𝑓𝑓 𝑦𝑦𝑓𝑓 𝑧𝑧𝑓𝑓 from chromosome 𝑥𝑥⃗. To guarantee a clinically acceptable reduction alignment, 𝑓𝑓 must be within a user-defined error boundary(≤ 0.5 𝑚𝑚𝑚𝑚 ). With this approach, it is possible to control the precision for specific regions on the bone. It allows not only a precise control of the reduction, but also the introduction of small deviations to the reconstruction target. Such deviations can support finding a better overall solution in case of contradicting targets by accepting small errors in the reduction alignment in favor of improving other objectives. For example, in case that the position of the fixation plate with respect to the 𝐴𝐴 -transformed bone fragment might not be 𝑓𝑓 within the acceptable solution range, the algorithm could allow a larger transformation error, e.g., 0.5 mm, and this small variation could already generate solutions with better fitting positions of the fixation plate. Control of osteotomy plane The parts of the bones that deviate from the anatomical target and that entail a surface gap among the fragments, are often referred to as bone protrusion (Figure 7A). The optimization of bone protrusion is very challenging because it is dependent on both, the osteotomy plane and the reduction alignment. We have decided to control the osteotomy cut position and direction implicitly by minimizing the bone protrusion among the generated fragments, in order to have a measurement that can be easily described to surgeons. As a positive side effect, the explicit optimization of the osteotomy plane is avoided, which has been proven to be challenging (Athwal et al., 2003; Schkommodau et al., 2005). The bone protrusion function 𝐵𝐵 𝑃𝑃 is defined as ( ) ( ) ( ) ( ) 𝐵𝐵 𝑃𝑃 𝐹𝐹 ,𝑅𝑅 = � 𝑊𝑊 𝑙𝑙𝑝𝑝 ,𝑝𝑝 − 𝛿𝛿 � 𝑙𝑙𝑝𝑝 ,𝑊𝑊𝑝𝑝 ,𝑊𝑊 𝑙𝑙𝑝𝑝 ,𝑅𝑅 � , (𝐄𝐄𝐄𝐄 .𝟕𝟕 ) 𝑚𝑚 𝑚𝑚 |𝐹𝐹 | 𝑝𝑝 ∈𝐹𝐹 𝑖𝑖 where 𝐹𝐹 is the point set of a bone fragment and 𝑅𝑅 is the reconstruction target. Windows thickness 𝑡𝑡 is set to � � � � ⃗ � � � � ⃗ 𝑡𝑡 = 4 𝑚𝑚𝑚𝑚 and 𝑡𝑡 and 𝑡𝑡 large enough to enclose 𝐹𝐹 and 𝑅𝑅 along 𝐶𝐶𝐶𝐶 and 𝐶𝐶𝐶𝐶 . 𝑦𝑦 𝑚𝑚 𝑧𝑧 𝑚𝑚 𝑧𝑧 In the optimization process, fitness function 𝑓𝑓 is defined as 𝑝𝑝 𝑟𝑟𝑝𝑝𝑚𝑚 𝑑𝑑 𝑚𝑚 𝐵𝐵 𝑃𝑃 (𝑃𝑃 ,𝑅𝑅 ) + 𝐵𝐵 𝑃𝑃 � 𝑃𝑃 𝐴𝐴 ,𝑅𝑅 � 𝑓𝑓 𝑓𝑓 = . (𝐄𝐄𝐄𝐄 .𝟖𝟖 ) 𝒑𝒑𝒑𝒑 An illustration of the bone protrusion calculation for 𝑷𝑷 is shown in Figure 7B-C. 𝒍𝒍𝒙𝒙 𝑑𝑑𝑑𝑑 Figure 7: Example of bone protrusion calculation. (A) Simplified drawing of the geometrical representation of bone protrusion after osteotomy cut and realignment. Gray cylinder represents the proximal fragment, orange cylinder represents 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 the reduced distal fragment. (B) Graphical explanation of the window function for fragment 𝑷𝑷 ; window function 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 ( ) 𝑾𝑾 𝒍𝒍𝒑𝒑 ,𝑷𝑷 is shown in blue, osteotomy plane is shown in cyan, the proximal fragment is shown in white and the 𝒙𝒙𝒑𝒑𝒍𝒍𝒑𝒑 ( ) reconstruction target is shown in transparent green. (C) Calculation of the bone protrusion 𝑷𝑷𝑩𝑩 𝑷𝑷 ,𝑹𝑹 , which can be 𝐩𝐩 interpreted as the calculation of the RMSE between the point sets 𝐏𝐏 (white points) and 𝑹𝑹 (green points) subject to the function window 𝑾𝑾 (see Eq. 7). Control of plate position Stability of the bone reduction and, consequently, successful healing, depends on an appropriate implantation of the fixation plate. In general, a minimal distance between fixation plate FP and bone surface is desirable. However, in orthopedic surgeries, the position of the fixation plate is often constrained by the approach site of the surgery (Klausmeyer and Mudgal, 2014) and by the intrinsic anatomy of the bones (muscle and ligament attachments and joints). To integrate this information into our algorithm, we have defined clinically feasible plate location regions according to all standard approach sites (Klausmeyer and Mudgal, 2014). These regions were annotated by clinical experts on an average mean model, generated by a statistical shape model (SSM) of the forearm (Mauler et al., 2017; Sepasian et al., 2014) (Figure 8A). The annotation was performed using the Scalismo mesh painting tool (Graphics and Vision Research Group, University of Basel, Switzerland). Region transfer from the mean model to each patient-specific model was achieved through a model fitting registration 𝑚𝑚 algorithm described by (Lüthi et al., 2018). Let 𝐵𝐵 𝑅𝑅 ∈ 𝑃𝑃 the i-th feasible bone region fitted to the pathological ( ) bone that can be retrieved by a function 𝐵𝐵 𝑅𝑅 𝑃𝑃 ,𝑒𝑒 . 𝑚𝑚 As a first pre-processing step, the fixation plate is coarsely aligned to each 𝐵𝐵 𝑅𝑅 of 𝑃𝑃 through ICP methods (Besl and McKay, 1992) and used as input for the optimization algorithm. In a second preprocessing step, an that should be brought into contact automatic identification is performed to determine plate model points 𝐹𝐹𝑃𝑃 𝑃𝑃𝐹𝐹 � � � � ⃗ with the bone surface. To this end, an inherent coordinate system 𝐶𝐶𝐶𝐶 for the fixation plate is calculated using a principal component analysis (PCA; (Jolliffe, 2011)) of FP. The eigenvectors of the PCA correspond to the 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 � � � � ⃗ Euclidian coordinates of the 𝐶𝐶𝐶𝐶 and the origin 𝑐𝑐𝑒𝑒 is given by the mean of 𝐹𝐹𝑃𝑃 . Without loss of generality, let 𝑝𝑝 𝑃𝑃𝐹𝐹 � � � � ⃗ the positive direction of the 𝐶𝐶𝐶𝐶 axis point towards the undersurface of the plate. A point of 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 is 𝑧𝑧 𝑚𝑚 ∗ 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 𝑃𝑃𝐹𝐹 � � � � ⃗ included in point set 𝐹𝐹𝑃𝑃 if (𝑝𝑝 − 𝑐𝑐𝑒𝑒 )∙ 𝐶𝐶𝐶𝐶 > 𝜀𝜀 𝑡𝑡 , where 𝑡𝑡 is the thickness of the fixation plate along 𝑚𝑚 𝑝𝑝 𝑧𝑧 𝑧𝑧 𝑧𝑧 𝑃𝑃𝐹𝐹 ∗ � � � � ⃗ 𝐶𝐶𝐶𝐶 . We empirically set 𝜀𝜀 = 0.85 , which worked for all plate types. An example of the identified 𝐹𝐹𝑃𝑃 point 𝑧𝑧 set is shown in Figure 8B. 𝐩𝐩𝐩𝐩𝐩𝐩 Figure 8: Feasible plate regions. (A) Feasible plate regions 𝑩𝑩 𝑹𝑹 (light yellow), encoded in the mean model of a radius 𝒎𝒎 SSM. The remaining colored regions are shown as indication of non-feasible plate regions (B) Surface points 𝑭𝑭𝑷𝑷 (in green) 𝑭𝑭𝑷𝑷 � � � � � ⃗ of the edge of the plate that should be in contact with the bone after registration. Coordinate axis 𝑪𝑪𝑪𝑪 is shown in black and 𝒛𝒛 is oriented towards the undersurface of the plate. 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 On each optimization step, the average distance between the fixation plate and the proximal fragment 𝑃𝑃 is calculated as 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 ( ) 𝐷𝐷 = �� 𝛿𝛿 � 𝐴𝐴 𝑝𝑝 ,𝑅𝑅𝐵𝐵 𝑃𝑃 ,𝑒𝑒 � − 𝑝𝑝 � , 𝑒𝑒 ℎ 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 (Eq. 9) ∗ 𝑝𝑝 𝑗𝑗 𝑗𝑗 𝑗𝑗 � 𝐹𝐹𝑃𝑃 � 𝑚𝑚 =1,..,|{𝐵𝐵 𝑅𝑅 }| 𝑖𝑖 𝑗𝑗 =1,…,� 𝑃𝑃𝐹𝐹 � ∗ ∗ where |𝐹𝐹𝑃𝑃 | is the cardinality of 𝐹𝐹𝑃𝑃 , 𝐴𝐴 is the 4x4 homogenous transformation matrix of the plate obtained 𝑝𝑝 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 𝑝𝑝 𝑟𝑟 𝑚𝑚𝑝𝑝 from chromosome 𝑥𝑥⃗, 𝐵𝐵 𝑅𝑅 (𝑃𝑃 ,𝑒𝑒 ) are the points of feasible region 𝑒𝑒 contained in 𝑃𝑃 , and 𝛿𝛿 is the already defined nearest neighbor function of Eq. 3. Similarly, the distance between the fixation plate 𝐹𝐹𝑃𝑃 and the distal 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 fragment 𝑃𝑃 is obtained by 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝐷𝐷 = �� 𝛿𝛿 � 𝐴𝐴 𝑝𝑝 ,𝐴𝐴 𝑅𝑅𝐵𝐵 (𝑃𝑃 ,𝑒𝑒 )� − 𝑝𝑝 � , where 𝑝𝑝 ∈ 𝐹𝐹𝑃𝑃 (Eq. 10) ∗ 𝑝𝑝 𝑗𝑗 𝑓𝑓 𝑗𝑗 𝑗𝑗 � 𝐹𝐹𝑃𝑃 � |{ }| 𝑚𝑚 =1,.., 𝐵𝐵 𝑅𝑅 𝑖𝑖 𝑗𝑗 =1,…,� 𝑃𝑃𝐹𝐹 � where 𝐴𝐴 is the 4x4 homogenous transformation matrix of the distal fragment obtained from chromosome 𝑥𝑥⃗. The 𝑓𝑓 average of both distance values is used for controlling the position and orientation of the plate, encoded into fitness function 𝑓𝑓 as 𝑝𝑝 𝑝𝑝𝑟𝑟 𝑚𝑚 𝑑𝑑 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑓𝑓 = � 𝐷𝐷 +𝐷𝐷 � ⁄2 (Eq. 11) Control of screw position and orientation The bone density of the patient bone considerably influences screw stability, as cortical and trabecular bone distribution differs along the bone (Kubiak et al., 2006; Reitman et al., 2004). The lengths of the screws play another important factor, as too short or too long screws might cause complications or injuries (Blecha et al., 2005; Spahn, 2004; Wall et al., 2012). Also, the position and direction of screws with respect to the bone density considerably affect the stability of the fixation (Dougherty et al., 2008; Weninger et al., 2010). Therefore, the bone density distribution will be considered for defining the screw positions. We have developed a method for generating a patient-specific bone density mask based on the original CT image 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 𝐼𝐼 containing the patient-specific normalized Hounsfield values (𝐼𝐼 ) (Figure 9A-B). Let Φ (𝑝𝑝 ) be a 𝐼𝐼 𝑚𝑚 transformation function transforming the 3D coordinates of a point 𝑝𝑝 to image space coordinates � ,𝐼𝐼 ,𝐼𝐼𝐼𝐼 � with 𝑚𝑚 𝑚𝑚 𝑦𝑦 𝑧𝑧 respect to an image 𝐼𝐼 . The outer cortical surface of the pathological bone model 𝑷𝑷 (orange overlay in Figure 9B) is 𝑃𝑃𝑇𝑇𝐿𝐿 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 rasterized into an empty image 𝐼𝐼 of same size as 𝐼𝐼 using the method described by (Pineda, 1988). A grey 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 value of 1000 𝑚𝑚𝑎𝑎𝑥𝑥 (𝐼𝐼 (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 )) is assigned to each rasterized voxel. Next, a dilation of size 1 pixel is (𝑚𝑚 ,𝑦𝑦 ,𝑧𝑧 ) 𝑃𝑃𝑇𝑇𝐿𝐿 applied to 𝐼𝐼 using a standard dilation kernel and the weighted bone-density mask shown in Figure 9C is defined 𝑒𝑒𝑒𝑒𝑒𝑒 𝐿𝐿 𝑊𝑊 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 ( ) ( ) ( ) 𝑀𝑀 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 = max� 𝑀𝑀 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ,𝐼𝐼 𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 � , (Eq. 12) 𝑃𝑃𝑇𝑇𝐿𝐿 𝑂𝑂 𝑟𝑟𝑚𝑚 𝑂𝑂 which returns the values of 𝐼𝐼 for (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ) lying on the cortical layer, values of 𝑀𝑀 (𝑥𝑥 ,𝑦𝑦 ,𝑧𝑧 ) lying inside the pathological bone (trabecular bone) and 0 otherwise. The masking process was performed in Matlab (MATLAB, 2017) using the Iso2Mesh toolbox (Fang and Boas, 2009). In the objective evaluation of the optimization process, penetration points 𝑒𝑒𝑛𝑛 and 𝑜𝑜𝑜𝑜 𝑡𝑡 (Figure 9D-E, cyan 𝑑𝑑 𝑑𝑑 𝑖𝑖 𝑖𝑖 points) are calculated between each screw 𝑒𝑒 and the bone fragments in reduced position. The penetration points 𝑚𝑚 � ⃗ are found by casting a ray formed by screw center A t (Eq. 4, with 𝑅𝑅 = I ) and screw axis h and p s 𝑑𝑑 4x4 s i 𝑖𝑖 i calculating the intersection with the bone as described in (Mount and Arya, 1998; Schneider and Eberly, 2002) . Afterwards, a sampling between 𝑒𝑒𝑛𝑛 and 𝑜𝑜𝑜𝑜 𝑡𝑡 is performed in 3D space along the screw axis using line equation 𝑑𝑑 𝑑𝑑 𝑖𝑖 𝑖𝑖 𝐿𝐿 � � =λ 𝑒𝑒𝑛𝑛 +λ � 𝑡𝑡 −𝑜𝑜𝑜𝑜 𝑒𝑒𝑛𝑛 � | λ = {0,0.1, … ,1} (Figure 9D, cyan and red spheres). The sampling 𝑃𝑃 𝑃𝑃 𝑑𝑑 𝑃𝑃 𝑑𝑑 𝑑𝑑 𝑃𝑃 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 𝑖𝑖 approach in 3D space is computationally more efficient than evaluating the entire density in the image space. The average weighted density value for each screw 𝑒𝑒 is obtained by 𝑚𝑚 𝑊𝑊 𝐷𝐷 𝑉𝑉 = � 𝑀𝑀 � Φ 𝑊𝑊 � 𝑙𝑙 � λ � � . (E� q. 13) 𝑃𝑃 𝐶𝐶 𝐶𝐶 𝑀𝑀 𝑖𝑖 𝑒𝑒 𝑒𝑒 � λ � 𝐶𝐶 𝑒𝑒 λ ={0,0.1,…,1} 𝐶𝐶 𝑒𝑒 The position and direction of the screws is controlled by the average of all the weighted density values 𝑁𝑁 𝑆𝑆 𝑖𝑖 𝑓𝑓 = � 𝐷𝐷 𝑉𝑉 (Eq. 14) 4 𝑃𝑃 𝑖𝑖 𝑁𝑁 𝑃𝑃 𝑖𝑖 𝑚𝑚 =1 𝐶𝐶𝑇𝑇 Figure 9: Fitness calculation for screw position and direction using patient-specific bone density information. (A) 3D 𝑶𝑶𝒑𝒑𝒎𝒎𝑶𝑶 models of the pathological radius (in orange) and ulna (in white); grey plane indicates a CT slice. (B) Slice of 𝑰𝑰 in coronal view. The radius is indicated by the semi-transparent orange mask and the outer cortical layer by the orange line. (C) 𝑾𝑾 ( ) Resulting weighted bone-density mask 𝑴𝑴 𝒙𝒙 ,𝒚𝒚 ,𝒛𝒛 . (D) Sampling calculation of insertion points for Screw 𝑪𝑪 ; the cyan 𝒎𝒎 spheres represent 𝒎𝒎 𝒎𝒎 and 𝒍𝒍 𝒐𝒐𝒕𝒕 and the red spheres represent the sampled points along the screw axis. (E) Coronal view of 𝒅𝒅 𝒅𝒅 𝒎𝒎 𝒎𝒎 𝑾𝑾 𝑴𝑴 (𝒙𝒙 ,𝒚𝒚 ,𝒛𝒛 ) showing the corresponding points of 𝑳𝑳 (𝛌𝛌 ) in image space. 𝑪𝑪 𝒎𝒎 2.3.3. Core Algorithm The core optimization algorithm is a modified version of the multi-stage multi-objective genetic optimization approach previously described in our work (Carrillo et al., 2017), capable of handling multiple contradicting objectives and nonlinear constraint. Once each planning objective has been parameterized using the real-valued chromosome 𝑥𝑥⃗, an optimal Pareto set 𝑌𝑌 of planning solutions can be found by 𝑚𝑚 𝑚𝑚 𝑚𝑚 { ∶ 𝑦𝑦 = min ( 𝑓𝑓 ( ) ( ) ( ) ) } : ℝ → ℝ , (Eq. 15) 𝑌𝑌 = 𝑦𝑦 𝜖𝜖 ℝ 𝑥𝑥⃗ , 𝑓𝑓 𝑥𝑥⃗ , … , 𝑓𝑓 𝑥𝑥⃗ | 𝑓𝑓 𝑚𝑚 ∈ Χ 1 2 𝑚𝑚 𝑚𝑚 𝑚𝑚 𝑚𝑚 where ℝ is the solution space generated by the fitness function 𝑓𝑓 (𝑥𝑥⃗), ℝ being the parameter space and Χ the 𝑚𝑚 𝑚𝑚 |𝑚𝑚⃗| 18 set of optimal chromosomes generated by the algorithm. In our case, ℝ = ℝ ≥ ℝ depending on the 𝑚𝑚 4 number of screws for each plate, and ℝ = ℝ where 𝑚𝑚 is the number of fitness functions. However, standard NSGA-II (Deb et al., 2002) is only able to find symmetrically optimized solutions with 𝑚𝑚 equally weighted optimization objectives on the solution space ℝ . In our optimization problem, each objective must have a different importance, according to its clinical relevance. Based on (Carrillo et al., 2017), we have applied the following weighting schema: ⁄ 𝑖𝑖 −k 𝑚𝑚⃗ 𝑃𝑃 𝑖𝑖 𝑒𝑒 (𝑥𝑥⃗) = ∑ 𝑒𝑒 . (Eq. 16) 𝑚𝑚 The weighting function 𝑒𝑒 increases the sparsity of solutions 𝑌𝑌 around the utopia point (Miettinen, 1999). The sparsity of 𝑓𝑓 (𝑥𝑥⃗) along the Pareto set is controlled by the constant 𝑘𝑘 , which represent the complementary 𝑚𝑚 𝑚𝑚 percentage of the desired sparsity. For instance, if 𝑘𝑘 = 30, 100− 𝑘𝑘 = 70% of solutions are to be found closer 𝑚𝑚 𝑚𝑚 to the utopia point. We have defined, together with surgeons, the following optimal weighting schema: Reduction alignment as the objective with the highest priority, followed by the osteotomy plane, the plate position and the screw purchase, respectively. An initialization and automatic verification of constraints and boundaries is done previous to each optimization stage. The distribution of the optimization process into different stages allows handling contradictive objectives by giving different importance on each stage to different objectives. Initialization Stage 1. An initial population 𝑋𝑋 , where the superscript denotes the optimization stage, of 200 chromosomes is randomly generated among the parameter range of chromosome 𝑥𝑥⃗ (Eq. 5). At each stage, we optimize over all parameters (all objectives) but considering only a subset of fitness functions. A summary of the algorithm parameters is given in Table 2. Stage 1. In this stage, the reduction alignment (𝑓𝑓 ), the osteotomy plane (𝑓𝑓 ), and the position of fixation plate 1 2 (𝑓𝑓 ) are optimized. Our algorithm is run with 3 objectives and weighted towards 𝑓𝑓 (𝑘𝑘 = 20, Eq. 16) and 𝑓𝑓 3 1 1 2 (𝑘𝑘 = 40, Eq. 16). This strategy permits to obtain solutions for the reduction alignment and bone protrusion close to the utopia point, but allowing a larger freedom on the plate position (𝑘𝑘 is set to 100, i.e., sparsity is unaffected). By forcing the algorithm to remain in the solution space where the reduction alignment (𝑓𝑓 ) and the osteotomy plane (𝑓𝑓 ) are close to the utopia point, contradictive solution that might engender the quality of 𝑓𝑓 or 2 1 𝑓𝑓 , are automatically neglected by the evolutionary process of the GA. The resulting Pareto set is stored in matrix 1 1∗ 𝑌𝑌 containing the best 70 solutions corresponding to the Pareto front. Accordingly, 𝑋𝑋 represents the corresponding set of parameters yielding 𝑌𝑌 . Table 2: Summary of algorithm parameters for each stage Parameter Stage 1 Stage 2 Input Population 𝑋𝑋 𝑋𝑋 Max # Generations 200 200 Population Size 200 200 Number of objectives 3 3 𝑓𝑓 , 𝑓𝑓 , 𝑓𝑓 𝑓𝑓 , 𝑓𝑓 , 𝑓𝑓 Fitness Functions 𝑓𝑓 (𝑥𝑥⃗) 1 2 3 1 3 4 𝑚𝑚 𝐾𝐾 = 20, 𝑘𝑘 = 40 𝐾𝐾 = 40 Weighting Schema(𝐾𝐾 ) 1 2 3 𝑚𝑚 1 1∗ 2 2∗ Output 𝑌𝑌 , 𝑋𝑋 𝑌𝑌 , 𝑋𝑋 Initialization Stage 2. The best 70 solutions from stage 1 are used for the initialization of the first generation of stage 2. To this end, the parameter range 𝑋𝑋 of the initial population of stage 2 are randomly initialized per component, such that each component value is constrained to lie between the maximum and minimum values 1∗ 1∗ 1∗ obtained from 𝑋𝑋 . Thus, 𝑥𝑥⃗ ∈ � 𝑚𝑚𝑒𝑒𝑛𝑛 �� ,𝑋𝑋𝑚𝑚𝑎𝑎𝑥𝑥 �� 𝑋𝑋 wher� e 𝑒𝑒 is the component i-th of the chromosome 𝑥𝑥⃗. 𝑚𝑚 𝑚𝑚 𝑚𝑚 The algorithm parameters of stage 2 are given in Table 2. Stage 2. In this stage, we optimize the reduction alignment (𝑓𝑓 ), the position of fixation plate (𝑓𝑓 ) and the 1 3 average screw purchase (𝑓𝑓 ). Our algorithm is weighted towards 𝑓𝑓 (𝑘𝑘 = 40, 𝑬𝑬𝑬𝑬 .𝟏𝟏𝟏𝟏) , to guarantee solutions 4 3 3 with an optimal fixation plate alignment, and to avoid deteriorating the reduction alignment. The remaining fitness functions are not weighted (𝑘𝑘 = 𝑘𝑘 = 100). The resulting Pareto set is stored in matrix 𝑌𝑌 , containing 1 2 2 the best 70 optimal chromosomes corresponding to the Pareto front. Thanks to our weighing schema, 𝑌𝑌 contains 2∗ 2 optimal solutions that are all within an acceptable clinical errors. 𝑋𝑋 was ranked according to the solutions 𝑌𝑌 2∗ 𝑃𝑃𝑝𝑝𝑙𝑙 by the best combined fitness among the 3 objectives of stage 2, yielding 𝑋𝑋𝑅𝑅 . The final output 𝑥𝑥⃗ of the algorithm is obtained by taking the solution with the best plate alignment among the ranked list � 𝑓𝑓 (𝑥𝑥� � ⃗)− 𝑚𝑚𝑒𝑒𝑛𝑛 (𝑓𝑓 (𝑥𝑥� � ⃗))� 3 3 𝑃𝑃𝑝𝑝𝑙𝑙 𝑥𝑥⃗ = min . 𝑥𝑥� � ⃗ ∈ 𝑅𝑅𝑋𝑋 (� � ⃗) (� � ⃗) � 𝑚𝑚𝑎𝑎𝑥𝑥 (𝑓𝑓 𝑥𝑥 − 𝑚𝑚𝑒𝑒𝑛𝑛 (𝑓𝑓 𝑥𝑥 � 3 3 3. Results 3.1. Datasets We have performed a qualitative validation (section 3.2) and a quantitative evaluation (section 3.3) on retrospective cases of malunited radii, which have been included in a large clinical trial about CA corrective osteotomy. From these consecutive cases, 36 cases were eligible according to the inclusion criteria given in Table 3. All 36 patients were treated at our orthopedics department between 2015 and 2017 and underwent navigated forearm osteotomy surgery through 3D preoperative planning and patient-specific instrumentation. Informed consent to use the patient data in an anonymized form for computer analyses was obtained, with ethical approvals KEK-ZH-Nr. 2015-0186 and BASEC-Nr. 2018-01608. The corresponding CT scans used for the generation of the 3D models were obtained according to a standard scanning protocol (slice thickness 1 mm; 120 kV; Philips Brilliance 64 CT, Philips Healthcare, The Netherlands). 3D preoperative planning of all patients were prospectively created in a manual fashion using a preoperative planning software (CASPA, Balgrist CARD AG, Switzerland) by the responsible hand surgeon, together with an expert planning engineer. These solutions were considered as the Gold Standard (GS). The baseline of the selected data is presented in Table 3. Table 3: Inclusion and exclusion criteria for clinical and technical evaluation of the clinical cases Inclusion Exclusion  Signed informed consent for data use  Incomplete CT data  Revision surgery  Age ≥ 16 years old  Radius osteotomy  Intra-articular corrective osteotomy  Availability of both, pathological and  Multiple osteotomies on one bone contralateral CT data  Complete 3D preoperative planning including fixation plate Table 4: Baseline data of the 36 cases included in the study Mean SD Age (years)* 33 ±14 Male Female Gender 11 25 Right Left Affected Side 19 17 Distal Radius Radius Shaft Location Malunion 14 22 Single Closing Opening Cut Wedge Wedge Osteotomy Type 11 8 17 Fixation Plate 9 different types * years at time of surgery 3.2. Qualitative Validation A preoperative planning solution (optimization solution; OA) was generated by our algorithm for each of the 36 consecutive cases in an automated fashion. The anatomical site and osteosynthesis implant were defined to be the same as in the GS solution for each case. A clinical survey was performed with 6 readers (2 senior hand surgeons, 1 board-qualified orthopedic surgeon, 1 resident surgeon and 2 expert planning engineers). In the survey, we asked each reader the question: “which of the two presented solutions would you implement in the surgery without further modification?”. The two solutions (GS and OA) were presented to the readers in a blinded fashion and in random order. The assessment for each case is shown in Figure 10A. OA solutions were chosen as the better solution in 55% of the times (Figure 10B). Out of the 55%, surgeons’ votes represented 38% and engineers’ votes 17%. Similarly, for the 45% of the GS votes, 29% account for the votes of the surgeons and the rest 16% represent the votes of the engineers. Figure 10: Results of the qualitative validation. (A) Voting of each reader and case. Each case was evaluated by 4 surgeons (readers 1-4) and 2 experienced planning engineers (readers 5-6). (B) Percentage distribution among the different voters for GS and OA solutions. In Figure 11 we present two examples where the OA solutions outperformed the GS solutions. In both cases, the reduction alignment and the position of the fixation plate obtained by the OA were clearly better than the GS solutions. Moreover, in both cases, the distance of the fixation plate with respect to the bone fragments has been decreased and undesired bone intersections have been solved by the OA. In these two cases, 100% of the readers favored the OA solutions. There were also cases where the GS solutions were preferred by the majority of the readers. Figure 12 shows two example cases where more than 50% of the readers favored the GS solutions. In the case shown in Figure 12A, surgeons preferred the GS solution due to a smaller bone protrusion with respect to the corresponding OA solution. In the case shown in Figure 12B, the GS solution has a better plate positioning, albeit a slightly worse reduction alignment compared to the OA solution. In 8 out of the 36 cases, the GS solution corresponded to a closing wedge osteotomy (C1, C5-C7, C10, C18, C20, and C21) because planners failed to manually find a clinically acceptable single-cut osteotomy. For all of these 8 cases, the optimization framework was rerun and the algorithm converged towards a single-cut solution by applying the corresponding constraints. A single-cut osteotomy is the most preferred clinical implementation of a corrective osteotomy. Figure 13 shows one of the OA single-cut solutions in comparison to the closing- wedge approach of the GS (Figure 13A). In the OA solution (Figure 13B), not only the osteotomy type was improved, but also the position of the fixation plate. Figure 11: Example cases for OA solutions that outperformed GS solutions. The proximal fragment of the pathological bone is shown in white and aligned to the reconstruction target, shown in green. The reduced distal fragments of the OA and GS solutions are shown in magenta and purple, respectively. (A) Case 21 depicts a radius shaft osteotomy where the OA solution was able to solve undesired plate-bone intersections. (B) Case 30 shows a distal radius osteotomy. The OA solution was able to find a better positioning of the fixation plate avoiding bone intersections and improving the bone-plate contact. Figure 12: Example cases where GS solutions were assessed to be better than OA solutions. The proximal fragment of the pathological bone is shown in white and aligned to the reconstruction target, shown in green. The reduced distal fragments of the OA and GS solutions are shown in magenta and purple, respectively. (A) Case 07 depicts a radius shaft osteotomy where the bone protrusion was assessed to be better in the GS solution. The bone protrusion in the OA solution was larger in comparison to the GS, albeit a slightly better reduction alignment (B) Case 26 shows a distal radius osteotomy where the position of the plate was assessed to be better in the GS solution. In the OA solution, the plate-bone distance was optimized for the proximal fragment, resulting in a small bone intersection on the distal fragment. Figure 13: Example of case 20, where the algorithm was able to find a single-cut osteotomy solution. The proximal fragment of the pathological bone is shown in white and the reconstruction target is shown in green. (A) Closing-wedge solution generated by the state-of-the-art planning (GS); realigned distal fragment is shown in magenta and the wedge to be removed to complete the osteotomy is shown in red. (B) Single-cut solution generated by our approach (OA); realignment of the distal fragment is shown in purple. Position of the fixation plate was also improved with respect to the GS solution. 3.3. Quantitative Evaluation Six different error measurements were evaluated for both, the GS and OA solutions, and for each of the 36 cases. Their average values, standard deviations, and ranges are reported in Table 5. (a) Reduction error, measured by Eq. 6, between the transformed distal fragment and the reconstruction target. A reduction error <= 0.5 mm is desired clinically. (b) Average distance between the fixation plate and the proximal bone fragment measured by Eq. 9. (c) Average distance between the fixation plate and the distal bone fragment measured by Eq. 10. (d) Average screw purchase of the proximal screws of the fixation plate, measured by Eq. 14. A larger screw purchase represents a better positioning of the screw. (e) Average screw purchase of the distal screws of the fixation plate, measured by Eq. 14. A larger screw purchase represents a better positioning of the screw. (f) Minimum distance between the screws of the fixation plate and the osteotomy plane measured by 𝑚𝑚𝑒𝑒𝑛𝑛 � 𝑛𝑛� ⃗ . (𝐶𝐶 − 𝑙𝑙𝑝𝑝 )� � � 𝑛𝑛� ⃗ � . A distance ≥2 mm is preferred. 𝑃𝑃 𝑝𝑝𝑙𝑙 𝑚𝑚 𝑝𝑝𝑙𝑙 𝑖𝑖 Table 5: Mean, standard deviation and range for the 6 error measurements used in the technical evaluation of the optimization solutions compared to the gold standard. Gold Standard (GS) Optimization Algorithm (OA) Measurement Mean ± SD Range Mean ± SD Range Alignment Error (mm) 1.47 1.89 0.08 , 9.56 0.94 0.88 0.07 , 3.55 Distance Plate-Bone Proximal (mm) 0.83 0.67 -1.04 , 2.64 0.57 0.44 -0.80 , 1.59 Distance Plate-Bone Distal (mm) 0.88 0.66 -1.05 , 2.03 0.64 0.44 -0.88 , 1.62 Screw Purchase Proximal (unit) 1.4 1.02 -0.86 , 2.90 1.98 0.74 -0.22 , 3.34 Screw Purchase Distal (unit) 1.15 0.89 -0.22 , 3.26 1.32 0.73 -0.06 , 2.69 Minimum Distance 4.91 3.05 0.45 , 9.59 5.75 2.53 0.80 , 12.34 Screw-Osteotomy Plane (mm) We have also performed an evaluation with respect to each planning objective by assessing the error measures for both GS and OA solutions. Figure 14 shows a set of box plots expressed as the relative improvement of the measures that were achieved by the OA solution. The blue baseline represents the GS value. In average, the alignment error was improved by 20% with respect to the GS among all the evaluated cases. The distance between the plate and the proximal and distal bone fragments has been improved by 26% and 31%, respectively. The average screw purchase for the proximal and distal screws was improved by 106% and 107%, respectively. Finally, an improvement of 93% in the average minimum distance between the fixation screws and the osteotomy plane was achieved by the OA solutions. Figure 14: Results of quantitative evaluation. Box plot (3/2 interquartile range whiskers) of the relative improvement of OA solutions in relation to the GS (blue vertical bold line) for each of the 36 cases. The black bold line indicates the median and the red star indicates the average improvement for each error measurement. The black crosses indicate the outliers. 3.4. Runtime The presented approach was implemented in Matlab R2017b (MATLAB, 2017). The average runtime of the algorithm comprising pre-processing, automatic diagnosis, verification and initialization of constraints, two optimization stages and generation of output solution was 85.38 min (±33.15 min). 4. Discussion In the last decade, the use of computer-assisted 3D preoperative planning for orthopedic surgeries has increased significantly due to its higher precision and to the capability of treating more complex pathologies (Murase et al., 2008; Nagy et al., 2008). However, the great effort required for generating preoperative planning solutions with current state-of-the-art approaches poses a bottleneck in the treatment of corrective osteotomies. In this work, we have presented a computer framework for the fully automatic calculation of preoperative planning solutions of osteotomies. Our approach augments the skills of the surgeons by providing improved solutions with a reliable computational tool. To the best of our knowledge, this automated preoperative planning framework represents the first to have been developed to generate ready-to-use preoperative planning solutions for corrective osteotomies of the forearm In orthopedic surgery, good surgical outcome is tightly associated with careful preoperative planning and precise surgical execution (Miyake et al., 2012b; Vlachopoulos et al., 2015; Vroemen et al., 2012). In order to demonstrate the clinical validity of our framework, we conducted a clinical study which confirmed the capabilities of the algorithm for generating ready-to-use preoperative solutions. The study showed that solutions generated by our algorithm were preferred in 55% of the time. In the remaining 45%, the quantitative evaluation confirmed still a high quality of the OA solution, with an error to the GS solutions of less than 1 mm for the reduction alignment, and less than 1.5 mm for the bone-plate distance. One clear advantage of our framework is the ability of generating solutions in an automated fashion, deducting human calculation and planning times from the clinical setting, in contrast to state-of-the-art methods, which are associated with high effort and clinical costs (Fürnstahl et al., 2016). The calculation of 3D preoperative planning solutions in an automatic way would represent a time reduction of almost 35% of the human workload per case, decreasing therefore significantly the related costs. Our institution performs more than 250 corrective osteotomies each year using a commercial service, in which the 3D preoperative planning is performed in a manual fashion. If we consider average diagnosis and planning costs of USD 600 per case, our approach would save USD 150’000 of yearly clinical costs in our institution alone, while improving also the quality of planning and patient treatment. Another advantage is the capability of calculating solutions with different trade-offs among the optimization objectives. Our algorithm can drastically improve one of the objectives by slightly decreasing the quality of other objectives. Such evaluations are not possible to obtain in a systematic way using the manual state-of-the-art approach (Bauer et al., 2017b; Vlachopoulos et al., 2015). Moreover, our method allows calculating single-cut osteotomies for arbitrary bone deformities including angular and translational deformities. In the presented case series, 8 single-cut OA solutions could be calculated where the real surgery was implemented with a closing wedge osteotomy. In 75% of these cases, readers validated the generated OA single-cut osteotomy solutions as a better preoperative planning than the original closing-wedge osteotomy proposed by a surgeon. Concerning the surgical execution, the generation of complex solutions for preoperative planning comes also at the cost of an increase in the complexity of the intraoperative navigation. At our institution, we have developed a state-of-the-art technique for intraoperative guidance using PSI, which are 3D-printed individually for each case. These PSI allow a high level of precision for executing the preoperative plan in the surgery (Fürnstahl et al., 2016; Omori et al., 2014; Rosseels et al., 2018). An example of the generated PSI for intraoperative navigation of a distal radius osteotomy is shown in Figure 15. The accuracy of PSI navigation for corrective osteotomies has been previously evaluated, showing an improve in the postoperative outcome (Asada et al., 2014; Bauer et al., 2017a; Byrne et al., 2017; Caiti et al., 2018; Farshad et al., 2017; Gouin et al., 2014; Kunz et al., 2013; Michielsen et al., 2018; Omori et al., 2014; Roner et al., 2018; Rosseels et al., 2018; Schweizer et al., 2016; Vlachopoulos et al., 2015) Figure 15: Example of translation of preoperative plan to the intraoperative situation by means of 3D-printed PSI. (A) Preoperative plan of a distal radius osteotomy. (B) PSI design for the osteotomy cut. (C) PSI design for the implant navigation and bone reduction. (D) Intraoperative navigation of the osteotomy cut. (E) Intraoperative navigation of the bone reduction and implant positioning. The black arrows indicate the PSI design that corresponds to each intraoperative situation. The presented framework has some technical and clinical limitations. In the first place, the number of cases included in this study was not enough to provide statistical relevance to the quantitative and qualitative analysis performed. However, the observed trend shows a promising potential of the capabilities of the developed framework. The quality of the solutions generated by the optimization algorithm was very similar for every case to the ones of the Gold Standard. The decision towards one or the other was determined by very specific factors such as a better angle for the osteotomy cut or a better position for the fixation plate. Nevertheless, this is also the case when treating patients with the Gold Standard method as the surgeon has to decide which solution will be the one transferred to the intraoperative situation by means of PSI (Roner et al., 2018; Rosseels et al., 2018; Schweizer et al., 2013; Victor and Premanathan, 2013; Vlachopoulos et al., 2015; Wong et al., 2017) In order to study the accuracy of surgical outcomes using automatic optimization methods, a prospective study with a statistical significant number of cases needs to be performed including a control group. Nevertheless, before such a clinical study can be approved and performed, the performance and safety of the algorithm needed to be validated. We have successfully achieved that in this paper through two different steps: (1) a quantitative evaluation showing that the algorithm is able to generate the expected solutions, optimizing all objectives, and (2) qualitative evaluation by a survey among surgeons and engineers showing that the algorithm produces solutions that the operating surgeon would execute in the practice. Secondly, the optimization is not able to take into account clinical decisions made by the surgeons based on patient history, for example, an unusual reduction alignment due to soft tissue pathologies (Vlachopoulos et al., 2015), or a different fixation technique due to osteoporotic bones. Also, the presence of a deformity in the surrounding bones has a significant impact in the final reduction. In our algorithm, this influence was taken into account through the reconstruction target, however in future works, the influence of neighboring bones should be also taken into account as a part of the optimization framework. Furthermore, surgeons have intraoperative techniques to improve poorly fitting plates such as manual plate-bending or bone reshaping. Such actions cannot be foreseen by the algorithm. This might be also a reason why the GS solution was preferred in some cases even with a poor fitting of the fixation plate. Example of that can be seen for cases 7, 17, and 20, where the GS solution was preferred by 83% of the readers. In the future, patient-specific osteosynthesis and osteotomy implants such as the one described by Caiti et al. (2018) and Omori et al. (2014) could be integrated into the optimization framework. Our approach is an important first step to turn automatic osteotomy planning into clinical practice. Due to the modularity of the algorithm, the approach could likely be extended to handle further anatomies (e.g., hip, humerus, clavicle, wrist bones) (Bugeau and Pérez, 2008; Chahla et al., 2016; Reising et al., 2013; Tschannen et al., 2016; Vlachopoulos et al., 2018). The positioning of the fixation plate and the calculation of the screw purchase could directly be adapted without significant efforts, but the calculation of the osteotomy plane and bone reduction would require the adaptation of the fitness functions and constraints. Our algorithm does not include yet the possibility to support simultaneous osteotomies on multiple locations (Belei et al., 2007; Fürnstahl et al., 2016; Schkommodau et al., 2005), and it does not include either the capability to plan intra-articular osteotomies (Schweizer et al., 2013). However, we consider that the strategy could be adapted by including extra stages to the optimization pipeline, at the cost of increasing the calculation times due to the linear structure of the optimization stages. One possibility to avoid the increase of calculation times is to implement a parallel design for the optimization stages, also allowing the capability for a back- feedback mechanism of the entire optimization. The presented approach is 18% faster in comparison to our previous work (Carrillo et al., 2017). However, the computational effort is still very high. One of the reasons for long calculation times is the implementation of the algorithm in a non-compiled programming language. Another possible source for the long calculation time is the lack of parallelization. We expect that the implementation of the entire pipeline in a compiled language and the use of graphical processing units (GPUs) for a customized parallelization of the approach will speed up the optimization framework to interactive rates. 5. Conclusions and Further work The presented framework is able to generate clinically feasible preoperative planning solutions in an automatic fashion. The key idea of the approach is the use of a multi-staged, multi-objective optimization pipeline for the concurrent optimization of the planning objectives. We demonstrated that the clinical and technical quality of the solutions generated by the algorithm can be of the same quality as the solutions created by experienced surgeons. Moreover, in more than 50% of the cases, our algorithm was able to generate solutions that are hardly possible to be found within a reasonable time by the trial and error method of the Gold Standard. This shows the potential of the method to improve the quality of the preoperative planning by considering a larger range of feasible preoperative solutions. The algorithm does not intend to replace the expertise of surgeons in the process; our aim is to reduce unnecessary human work load. We seek to provide surgeons with automated tools for patient diagnosis and treatment, which can substantially decrease associated clinical times and costs, and generate better clinical outputs. Future works will focus on more flexible and less-invasive techniques for intraoperative guidance than patient-specific instrument (PSI) approaches such as augmented reality techniques in combination with learning-based methods for more comprehensive surgical guidance. Future solutions could also include a completely adaptive framework, were the optimal preoperative planning solution is adapted in real time according to the intraoperative situation. Conflict of interest None of the authors have any financial or personal relationships with other people or organizations that could bias their work. Acknowledgements This work has been funded through a Balgrist Foundation grant and a highly specialized medicine grant (HSM2) of the Canton of Zurich. We would also like to acknowledge the support and valuable input from the surgeons and planning engineers at the department of orthopaedics of Balgrist University Hospital. References Andress, S., Johnson, A., Unberath, M., Winkler, A.F., Yu, K., Fotouhi, J., Weidert, S., Osgood, G., Navab, N., 2018. On- the-fly augmented reality for orthopedic surgery using a multimodal fiducial. J Med Imaging (Bellingham) 5, 021209. Asada, S., Mori, S., Matsushita, T., Nakagawa, K., Tsukamoto, I., Akagi, M., 2014. Comparison of MRI- and CT-based patient-specific guides for total knee arthroplasty. The Knee 21, 1238-1243. Athwal, G.S., Ellis, R.E., Small, C.F., Pichora, D.R., 2003. Computer-assisted distal radius osteotomy1. The Journal of Hand Surgery 28, 951-958. Bauer, A.S., Storelli, D.A.R., Sibbel, S.E., Mccarroll, H.R., Lattanza, L.L., 2017a. Preoperative computer simulation and patient-specific guides are safe and effective to correct forearm deformity in children. Journal of Pediatric Orthopaedics 37, 504-510. Bauer, D.E., Zimmermann, S., Aichmair, A., Hingsammer, A., Schweizer, A., Nagy, L., Fürnstahl, P., 2017b. Conventional Versus Computer-Assisted Corrective Osteotomy of the Forearm: a Retrospective Analysis of 56 Consecutive Cases. The Journal of Hand Surgery 42, 447-455. Belei, P., Schkommodau, E., Frenkel, A., Mumme, T., Radermacher, K., 2007. Computer-assisted single- or double-cut oblique osteotomies for the correction of lower limb deformities. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 221, 787-800. Besl, P.J., McKay, N.D., 1992. Method for registration of 3-D shapes, Sensor Fusion IV: Control Paradigms and Data Structures. International Society for Optics and Photonics, pp. 586-607. Bilic, R., Zdravkovic, V., Boljevic, Z., 1994. Osteotomy for deformity of the radius. Computer-assisted three-dimensional modelling. Journal of Bone & Joint Surgery, British Volume 76-B, 150-154. Blecha, L.D., Zambelli, P.Y., Ramaniraka, N.A., Bourban, P.E., Månson, J.A., Pioletti, D.P., 2005. How plate positioning impacts the biomechanics of the open wedge tibial osteotomy; A finite element analysis. Computer Methods in Biomechanics and Biomedical Engineering 8, 307-313. Bugeau, A., Pérez, P., 2008. Track and Cut: Simultaneous Tracking and Segmentation of Multiple Objects with Graph Cuts. EURASIP Journal on Image and Video Processing 2008, 317278. Byrne, A.-M., Impelmans, B., Bertrand, V., Van Haver, A., Verstreken, F., 2017. Corrective Osteotomy for Malunited Diaphyseal Forearm Fractures Using Preoperative 3-Dimensional Planning and Patient-Specific Surgical Guides and Implants. Journal of Hand Surgery 42, 836.e831-836.e812. Caiti, G., Dobbe, J.G.G., Loenen, A.C.Y., Beerens, M., Strackee, S.D., Strijkers, G.J., Streekstra, G.J., 2018. Implementation of a semiautomatic method to design patient-specific instruments for corrective osteotomy of the radius. International Journal of Computer Assisted Radiology and Surgery. Carrillo, F., Vlachopoulos, L., Schweizer, A., Nagy, L., Snedeker, J., Fürnstahl, P., 2017. A Time Saver: Optimization Approach for the Fully Automatic 3D Planning of Forearm Osteotomies, Medical Image Computing and Computer-Assisted Intervention − MICCAI 2017. Springer International Publishing, Quebec, Canada, pp. 488-496. Chahla, J., Dean, C.S., Mitchell, J.J., Moatshe, G., Serra Cruz, R., LaPrade, R.F., 2016. Medial Opening Wedge Proximal Tibial Osteotomy. Arthroscopy Techniques 5, e919-e928. Criminisi, A., Shotton, J., 2013. Decision forests for computer vision and medical image analysis. Springer Science & Business Media. Deb, K., Agrawal, S., Pratap, A., Meyarivan, T., 2000. A fast elitist non-dominated sorting genetic algorithm for multi- objective optimization: NSGA-II. Lecture notes in computer science 1917, 849-858. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. Evolutionary Computation, IEEE Transactions on 6, 182-197. Dobbe, J.G.G., Strackee, S.D., Schreurs, A.W., Jonges, R., Carelsen, B., Vroemen, J.C., Grimbergen, C.A., Streekstra, G.J., 2011. Computer-Assisted Planning and Navigation for Corrective Distal Radius Osteotomy, Based on Pre- and Intraoperative Imaging. Biomedical Engineering, IEEE Transactions on 58, 182-190. Dobbe, J.G.G., Strackee, S.D., Streekstra, G.J., 2017. Minimizing the translation error in the application of an oblique single- cut rotation osteotomy: Where to cut? IEEE Transactions on Biomedical Engineering PP, 1-1. Dobbe, J.G.G., Vroemen, J.C., Strackee, S.D., Streekstra, G.J., 2014. Patient-specific distal radius locking plate for fixation and accurate 3D positioning in corrective osteotomy. Strategies in Trauma and Limb Reconstruction 9, 179-183. Dougherty, P.J., Kim, D.-G., Meisterling, S., Wybo, C., Yeni, Y., 2008. Biomechanical Comparison of Bicortical Versus Unicortical Screw Placement of Proximal Tibia Locking Plates: A Cadaveric Model. Journal of Orthopaedic Trauma 22, 399- Eck, M., Hoschek, J., Weber, U., 1990. Three-dimensional determination of an oblique osteotomy in the hip by mathematical optimization fulfilling some anatomical demands. Journal of biomechanics 23, 1061-1067. Esfandiari, H., Newell, R., Anglin, C., Street, J., Hodgson, A.J., 2018. A deep learning framework for segmentation and pose estimation of pedicle screw implants based on C-arm fluoroscopy. International Journal of Computer Assisted Radiology and Surgery 13, 1269-1282. Fang, Q., Boas, D.A., 2009. Tetrahedral mesh generation from volumetric binary and grayscale images, Biomedical Imaging: From Nano to Macro, 2009. ISBI'09. IEEE International Symposium on. Ieee, pp. 1142-1145. Farshad, M., Betz, M., Farshad-Amacker, N.A., Moser, M., 2017. Accuracy of patient-specific template-guided vs. free-hand fluoroscopically controlled pedicle screw placement in the thoracic and lumbar spine: a randomized cadaveric study. European Spine Journal 26, 738-749. Fernandez, D.L., 1982. Correction of post-traumatic wrist deformity in adults by osteotomy, bone-grafting, and internal fixation. J Bone Joint Surg Am 64, 1164-1178. Fürnstahl, P., 2010. Computer-assisted planning for orthopedic surgery. Diss., Eidgenössische Technische Hochschule ETH Zürich, Nr. 19102, 2010. Fürnstahl, P., Fuchs, T., Schweizer, A., Nagy, L., Székely, G., Harders, M., 2008. Automatic and robust forearm segmentation using graph cuts, Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on. IEEE, pp. 77-80. Fürnstahl, P., Schweizer, A., Graf, M., Vlachopoulos, L., Fucentese, S., Wirth, S., Nagy, L., Szekely, G., Goksel, O., 2016. Surgical Treatment of Long-Bone Deformities: 3D Preoperative Planning and Patient-Specific Instrumentation, In: Zheng, G., Li, S. (Eds.), Computational Radiology for Orthopaedic Interventions. Springer International Publishing, Cham, pp. 123- Gosse, F., Brack, C., Götte, H., Roth, M., Rühmann, O., Schweikard, A., Vahldiek, M., 1997. Roboterunterstützung in der KnieendoprothetikRobotic-assisted total knee replacement. Der Orthopäde 26, 258-266. Gouin, F., Paul, L., Odri, G.A., Cartiaux, O., 2014. Computer-Assisted Planning and Patient-Specific Instruments for Bone Tumor Resection within the Pelvis: A Series of 11 Patients. Sarcoma 2014, 1-9. Jolliffe, I., 2011. Principal component analysis, International encyclopedia of statistical science. Springer, pp. 1094-1096. Kawakami, H., Sugano, N., Nagaoka, T., Hagio, K., Yonenobu, K., Yoshikawa, H., Ochi, T., Hattori, A., Suzuki, N., 2002. 3D analysis of the alignment of the lower extremity in high tibial osteotomy, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, pp. 261-267. Klausmeyer, M.A., Mudgal, C., 2014. Exposure of the Forearm and Distal Radius. Hand Clinics 30, 427-433. Kubiak, E.N., Fulkerson, E., Strauss, E., Egol, K.A., 2006. The evolution of locked plates. JBJS 88, 189-200. Kunz, M., Ma, B., Rudan, J.F., Ellis, R.E., Pichora, D.R., 2013. Image-Guided Distal Radius Osteotomy Using Patient- Specific Instrument Guides. The Journal of Hand Surgery 38, 1618-1624. Lorensen, W.E., Cline, H.E., 1987. Marching cubes: A high resolution 3D surface construction algorithm, ACM siggraph computer graphics. ACM, pp. 163-169. Lüthi, M., Gerig, T., Jud, C., Vetter, T., 2018. Gaussian Process Morphable Models. IEEE Transactions on Pattern Analysis and Machine Intelligence 40, 1860-1873. MATLAB, 2017. 9.3.0.713579 (R2017b). The MathWorks Inc., Natick, Massachusetts. Mauler, F., Langguth, C., Schweizer, A., Vlachopoulos, L., Gass, T., Lüthi, M., Fürnstahl, P., 2017. Prediction of Normal Bone Anatomy for the Planning of Corrective Osteotomies of Malunited Forearm Bones Using a Three-Dimensional Statistical Shape Model. Journal of orthopaedic research: official publication of the Orthopaedic Research Society 35, 2630- Menetrey, J., Paul, M., 2004. Möglichkeiten der computergestützten Navigation bei kniegelenknahen Osteotomien. Der Orthopäde 33, 224-228. Merle d’Aubigné, R., Descamps, L., 1952. L’ostéotomie plane oblique dans la correction des deformations des membres. Bull Mem Arch Chirurgie 8, 271-276. Michielsen, M., Van Haver, A., Bertrand, V., Vanhees, M., Verstreken, F., 2018. Corrective osteotomy of distal radius malunions using three-dimensional computer simulation and patient-specific guides to achieve anatomic reduction. European journal of orthopaedic surgery & traumatology : orthopedie traumatologie 28, 1531-1535. Miettinen, K., 1999. Nonlinear Multiobjective Optimization. Springer US. Miyake, J., Moritomo, H., Kataoka, T., Murase, T., Sugamoto, K., 2012a. In vivo three-dimensional motion analysis of chronic radial head dislocations. Clinical Orthopaedics and Related Research® 470, 2746-2755. Miyake, J., Murase, T., Moritomo, H., Sugamoto, K., Yoshikawa, H., 2011. Distal Radius Osteotomy with Volar Locking Plates Based on Computer Simulation. Clinical Orthopaedics and Related Research® 469, 1766-1773. Miyake, J., Murase, T., Oka, K., Moritomo, H., Sugamoto, K., Yoshikawa, H., 2012b. Three-Dimensional Corrective Osteotomy for Malunited Diaphyseal Forearm Fractures Using Custom-Made Surgical Guides Based on Computer Simulation. JBJS Essential Surgical Techniques 2, e24. Mount, D.M., Arya, S., 1998. ANN: library for approximate nearest neighbour searching. Murase, T., Oka, K., Moritomo, H., Goto, A., Yoshikawa, H., Sugamoto, K., 2008. Three-Dimensional Corrective Osteotomy of Malunited Fractures of the Upper Extremity with Use of a Computer Simulation System. The Journal of Bone & Joint Surgery 90, 2375-2389. Nagy, L., Jankauskas, L., Dumont, C.E., 2008. Correction of forearm malunion guided by the preoperative complaint. Clinical orthopaedics and related research 466, 1419-1428. Omori, S., Murase, T., Kataoka, T., Kawanishi, Y., Oura, K., Miyake, J., Tanaka, H., Yoshikawa, H., 2014. Three‐ dimensional corrective osteotomy using a patient‐ specific osteotomy guide and bone plate based on a computer simulation system: accuracy analysis in a cadaver study. The International Journal of Medical Robotics and Computer Assisted Surgery 10, 196-202. Pineda, J., 1988. A parallel algorithm for polygon rasterization, ACM SIGGRAPH Computer Graphics. ACM, pp. 17-20. Reising, K., Strohm, P.C., Hauschild, O., Schmal, H., Khattab, M., Südkamp, N.P., Niemeyer, P., 2013. Computer-assisted navigation for the intraoperative assessment of lower limb alignment in high tibial osteotomy can avoid outliers compared with the conventional technique. Knee Surgery, Sports Traumatology, Arthroscopy 21, 181-188. Reitman, C., Nguyen, L., Fogel, G., 2004. Biomechanical evaluation of relationship of screw pullout strength, insertional torque, and bone mineral density in the cervical spine. Journal of spinal disorders & techniques 17, 306-311. Roner, S., Carrillo, F., Vlachopoulos, L., Schweizer, A., Nagy, L., Fuernstahl, P., 2018. Improving accuracy of opening- wedge osteotomies of distal radius using a patient-specific ramp-guide technique. BMC musculoskeletal disorders 19, 374. Roner, S., Vlachopoulos, L., Nagy, L., Schweizer, A., Fürnstahl, P., 2017. Accuracy and Early Clinical Outcome of 3- Dimensional Planned and Guided Single-Cut Osteotomies of Malunited Forearm Bones. The Journal of Hand Surgery 42, 1031.e1031-1031.e1038. Roscoe, L., 1988. Stereolithography interface specification. America-3D Systems Inc 27. Rosseels, W., Herteleer, M., Sermon, A., Nijs, S., Hoekstra, H., 2018. Corrective osteotomies using patient-specific 3D- printed guides: a critical appraisal. European Journal of Trauma and Emergency Surgery. Sangeorzan, B.P., Judd, R.P., Sangeorzan, B.J., 1989. Mathematical analysis of single-cut osteotomy for complex long bone deformity. Journal of Biomechanics 22, 1271-1278. Schenk, P., Vlachopoulos, L., Hingsammer, A., Fucentese, S.F., Fürnstahl, P., 2016. Is the contralateral tibia a reliable template for reconstruction: a three-dimensional anatomy cadaveric study. Knee Surg Sports Traumatol Arthrosc. Scheufler, K.-M., Franke, J., Eckardt, A., Dohmen, H., 2011. Accuracy of Image-Guided Pedicle Screw Placement Using Intraoperative Computed Tomography-Based Navigation With Automated Referencing, Part I: Cervicothoracic Spine. Neurosurgery 69, 782-795. Schkommodau, E., Frenkel, A., Belei, P., Recknagel, B., Wirtz, D.-C., Radermacher, K., 2005. Computer-assisted optimization of correction osteotomies on lower extremities. Computer Aided Surgery 10, 345-350. Schneider, P., Eberly, D.H., 2002. Geometric tools for computer graphics. Elsevier. Schweizer, A., Fürnstahl, P., Harders, M., Székely, G., Nagy, L., 2010. Complex radius shaft malunion: osteotomy with computer-assisted planning. Hand 5, 171-178. Schweizer, A., Fürnstahl, P., Nagy, L., 2013. Three-dimensional correction of distal radius intra-articular malunions using patient-specific drill guides. The Journal of hand surgery 38, 2339-2347. Schweizer, A., Mauler, F., Vlachopoulos, L., Nagy, L., Fürnstahl, P., 2016. Computer-assisted 3-dimensional reconstructions of scaphoid fractures and nonunions with and without the use of patient-specific guides: early clinical outcomes and postoperative assessments of reconstruction accuracy. The Journal of hand surgery 41, 59-69. Sepasian, N., Van de Giessen, M., Dobbe, I., Streekstra, G., 2014. Bone reposition planning for corrective surgery using statistical shape model: assessment of differential geometrical features, Bayesian and grAphical Models for Biomedical Imaging. Springer, pp. 49-60. Spahn, G., 2004. Complications in high tibial (medial opening wedge) osteotomy. Archives of Orthopaedic and Trauma Surgery 124, 649-653. Subburaj, K., Ravi, B., Agarwal, M., 2010. Computer-aided methods for assessing lower limb deformities in orthopaedic surgery planning. Computerized Medical Imaging and Graphics 34, 277-288. Tschannen, M., Vlachopoulos, L., Gerber, C., Székely, G., Fürnstahl, P., 2016. Regression forest-based automatic estimation of the articular margin plane for shoulder prosthesis planning. Medical Image Analysis 31, 88-97. Vatti, B.R., 1992. A generic solution to polygon clipping. Communications of the ACM 35, 56-63. Victor, J., Premanathan, A., 2013. Virtual 3D planning and patient specific surgical guides for osteotomies around the knee: a feasibility and proof-of-concept study. Bone Joint J 95-B, 153-158. Vlachopoulos, L., Carrillo, F., Gerber, C., Székely, G., Fürnstahl, P., 2017. A Novel Registration-Based Approach for 3D Assessment of Posttraumatic Distal Humeral Deformities. JBJS 99, e127. Vlachopoulos, L., Schweizer, A., Graf, M., Nagy, L., Fürnstahl, P., 2015. Three-dimensional postoperative accuracy of extra- articular forearm osteotomies using CT-scan based patient-specific surgical guides. BMC musculoskeletal disorders 16, 1. Vlachopoulos, L., Székely, G., Gerber, C., Fürnstahl, P., 2018. A scale-space curvature matching algorithm for the reconstruction of complex proximal humeral fractures. Medical Image Analysis 43, 142-156. Vroemen, J.C., Dobbe, J.G.G., Jonges, R., Strackee, S.D., Streekstra, G.J., 2012. Three-Dimensional Assessment of Bilateral Symmetry of the Radius and Ulna for Planning Corrective Surgeries. The Journal of Hand Surgery 37, 982-988. Wall, L.B., Brodt, M.D., Silva, M.J., Boyer, M.I., Calfee, R.P., 2012. The Effects of Screw Length on Stability of Simulated Osteoporotic Distal Radius Fractures Fixed With Volar Locking Plates. The Journal of Hand Surgery 37, 446-453. Weninger, P., Dall'Ara, E., Leixnering, M., Pezzei, C., Hertz, H., Drobetz, H., Redl, H., Zysset, P., 2010. Volar fixed-angle plating of extra-articular distal radius fractures—a biomechanical analysis comparing threaded screws and smooth pegs. Journal of Trauma and Acute Care Surgery 69, E46-E55. Wong, K., Paul, L., Gerbers, J., 2017. PATIENT SPECIFIC INSTRUMENTS CAN ACHIEVE A BETTER SURGICAL ACCURACY THAN NAVIGATION ASSISTANCE IN JOINT-PRESERVING SURGERY OF THE KNEE JOINT: A CADAVERIC COMPARATIVE STUDY. CAOS 1, 104-108. Zdravkovic, V., Bilic, R., 1990. Computer-assisted preoperative planning (CAPP) in orthopaedic surgery. Computer Methods and Programs in Biomedicine 32, 141-146. Vitae: Fabio Carrillo received his master degree in electronic engineering from Universidad Simon Bolivar in Caracas, Venezuela, with focus on automation. In 2015, he received his MSc degree in Robotics, Systems and Control from ETH Zürich, with a deep research interest in biomedical engineering. He joined the Computer Assisted Research and Development (CARD) Group at the University Hospital Balgrist in Zurich in September 2015, where he performs his doctoral studies, together with the Laboratory for Orthopaedic Biomechanics of the Department of Health Sciences and Technology at the ETH Zürich. Simon Roner is an orthopaedic surgeon resident working as a research associate in computer-assisted surgery. He graduated from the medical school at Zurich University in 2012. After completion of Swiss board exams in orthopaedic surgery, he will pursue a specialization in hand and upper limb surgery. Marco von Atzigen received his Master’s degree in Mechanical Engineering from ETH Zurich with focus on Robotics, Systems and Control in July 2017. He joined CARD in October 2017 as a software developer working on machine learning and augmented reality. Since December 2018, he is pursuing his doctoral studies in deep learning for medical imaging and scene understanding between CARD, the spine surgery team at the University Hospital Balgrist Zürich and the Laboratory for Orthopaedic Biomechanics of the Department of Health Sciences and Technology at the ETH Zürich. Andreas Schweizer graduated in Medicine from the University of Zurich. He received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 2003 and of Hand Surgery in 2005. He completed his MD thesis at Institute of Anatomy, University of Bern and his state doctorate at the University of Zurich. He is currently working as a consultant hand surgeon at the Balgrist University Hospital in Zurich, Switzerland and is specialized in 3D planning of orthopedic surgeries of upper extremity. Ladislav Nagy graduated in Medicine from the University of Zurich in 1982. He received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 1990 and of Hand Surgery in 1992. He is titular professor of the faculty of medicine of the University of Zurich since 2011. He is the current chief of the hand surgery of the university hospital Balgrist and work as a consultant specialist in 3D planning of orthopedic surgeries of the hand. Lazaros Vlachopoulos graduated in Medicine from the RWTH Aachen in Germany in 2004 and completed his MD thesis at the Institute of Human Genetics, RWTH Aachen. Afterwards, he completed his residency in orthopaedic surgery in Germany and Switzerland and received the Swiss Board Certification of Orthopaedic Surgery and Traumatology in 2012. In 2017 he received his PhD in medical image analysis from the ETH Zurich, Switzerland. He is currently working as a consultant orthopaedic surgeon at the Balgrist University Hospital in Zurich, Switzerland, specialized in computer-assisted orthopaedic surgery. Jess Snedeker is Associate Professor of Orthopaedic Biomechanics, holding a professorial chair at both the ETH Zurich (Department of Health Sciences and Technology) and the Medical Faculty of the University of Zurich (Department of Orthopaedics). He heads the division of experimental research at the University Hospital Balgrist, and also serves as the chief scientific officer of the Balgrist Campus, designated in 2017 by the Swiss Secretariat for Education, Research, and Innovation as a “Research Infrastructure of National Relevance”. Philipp Fürnstahl received the MSc degree in technical mathematics and information procession from the Technical University of Graz, Austria, in 2005. In 2010, he received the PhD degree in medical image analysis from the ETH Zurich, Switzerland, for his research in the field of computer-assisted preoperative surgery planning. Philipp Fürnstahl is currently head of the Computer Assisted Research and Development (CARD) Group at the University Hospital Balgrist in Zurich, Switzerland, where he is involved in the development of patient-specific solutions for orthopaedics surgeries.

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Sep 16, 2019

References