In this study, a computerized preoperative selection method is proposed for the treatment of femoral fractures with intramedullary nailing. The processing flow is as follows (Figure 1): First, the CT images of a femur are segmented and a 3D model is reconstructed. A centerline is extracted from the interested region of the MC. Second, the CT images of an IMN are segmented and its 3D model is reconstructed. The whole centerline of the IMN is extracted. Third, an entry point is preliminarily selected. The path from the entry point to the end of MC centerline is determined. Fourth, the centerline of the MC isthmus and the centerline of the IMN are registered. An optimal registration scheme is established to reduce the IMN surface points located outside the MC. Fifth, the IMN is registered to the MC based on the transformation obtained from the fourth step. Sixth, the operations of steps 2 to 5 are performed on each IMN. A suitable IMN can be determined with optimal orientation or alignment and entry point. Seventh, the distal nail hole can be locked accurately guided by a navigation system.

For both men and women, the left and right lower limbs maintain almost the same length, which is beneficial to maintain the balance of human bodies, and there is no significant difference between the left and right axes of the lower limbs under CT examination (Ferras-Tarrago et al., 2020). To facilitate the extraction of a femur centerline, images of the healthy side bone are used in this study. The healthy side bone and a candidate IMN are imaged by a CT scanner, and the DICOM sequence images are obtained. The Marching Cubes algorithm (Masala et al., 2013) is employed to reconstruct surface models of the femur (Figure 1A) and the IMN (Figure 1C) in three dimensions. The IMN is processed first. The centerline of the IMN is extracted by a homotopic thinning algorithm (Kerschnitzki et al., 2013). As a result, a central axis with a single voxel thickness is obtained, which is homotopic with the original image (Figure 2A). Then, the central axis is fitted as a curve line. A serial of control points with 1 mm interval are located on the curve line. For each control point, a fitting circle of IMN cross-section is placed. The diameter of the fitting circle and the curvature of IMN centerline around the control point are also recorded (Figure 2B).

Next, a 3D model of the healthy femur is imported into a 3D model processing software (Geomagic Control 2014). All the surface points of femur are removed manually, and the 3D model of the MC remains. Since the proximal end of femur includes the femoral head and the greater trochanter, and the distal end includes the condylar bone. The irregular shapes of both ends of the femur will affect the accuracy of extracting MC centerline. We only extract the middle tubular part as the interested part (Figure 2C). The centerline of the interested part is also extracted by the same homotopic thinning algorithm and control points are produced. The diameter of the fitting circle on the cross-section of the MC at each control point and the curvature of the MC centerline here are both recorded (Figure 2D).

Since the diameters of the distal part of an IMN are almost the same, while the MC is thick at both ends and thin in the middle. Besides, the curvature changes obviously in the middle narrow part of MC (marked as the MC isthmus) (Figure 2E). Therefore, the registration between the MC isthmus and the IMN is the key. Only the centerline of the interested part of MC is extracted in the previous step, so there is still a gap between the extracted centerline end and the apex of the greater trochanter. In order to locate an entry point and fill the gap, the intersection angle between the femoral neck’s axis and the shaft’s axis in 3D space, the femoral radius, the head radius, and the femoral length were measured. According to the method in Zhao et al. (2015) (Figure 1D), a possible area of entry points was determined and we manually select twelve candidate entry points. An interpolation method is used to place the control points between the candidate entry points and the starting point of interested part centerline respectively. Twelve central axes from the candidate entry points to the endpoint of interested part centerline which is near the intercondylar fossa are obtained. Next, we will find the most suitable one from the 12 central axes. If the control points on the one of central axes are used as the target point set and the control points on the IMN centerline as the source point set, a registration based on Iterative Closest Point (ICP) algorithm can be performed. Because ICP is easy to fall into the local optimal solution (Sobreira et al., 2019), this study combines Normal Distributions Transform (NDT) and ICP algorithm for registration (Shi et al., 2018). After registration, three influencing factors are calculated: curvature difference, diameter difference, and distance difference. The curvature difference between jth (j = 1, 2, … 12) central axis and the IMN centerline is calculated as Eq. 1: C j = ∑ i = 1 k C n a i l i − C c a v i t y i (1) where k is the number of control points of the IMN centerline, C n a i l i and C c a v i t y i are curvature at ith control points of the IMN centerline and the MC centerline. C Max and C Min denote the maximum value and minimum value of the set C 1 , C 2 , … C 12 respectively. The normalization C j M of C j is processed as follows: C j M = C j − C Min C Max − C Min (2)

The same operation in the calculation of C j M is performed on the diameter difference to obtain D j M . D j M represents the normalization of the diameter difference between jth central axis and the IMN centerline. The distance difference between jth central axis and the IMN centerline is calculated as Eq. 3: E j = ∑ i = 1 k E i (3) where E i is the Euclidean distance between ith control point of the IMN centerline and ith control point of the MC centerline. The normalization E j M of E j is calculated using the method in Eq. 2. The fitness of jth central axis registered with the IMN centerline is calculated as Eq. 4: φ j = C j M 3 + D j M 3 + E j M 3 (4)

The fitness φ varies between 0 and 1. The smaller the φ is, the more appropriate the central axis fits the IMN centerline. The most appropriate central axis is found out, and its corresponding entry point is used as the initial entry point.

In this part, a plane projection method is employed to calculate the surface points of the IMN laid outside the MC. According to the distribution of compenetration points, an iterative optimization of nail orientation, based on more and more extended volume regions including the isthmus is performed. Firstly, the narrowest part of the MC is calculated, and the initial length is set to 90 mm, which is used as the MC isthmus. Then, the centerline of the MC isthmus is intercepted. A centerline with the same length at the corresponding position of IMN is also intercepted. The two centerlines are registered by NDT and ICP algorithm, and a spatial transformation matrix T is obtained. Based on the transformation matrix T, the IMN is registered with the MC. We use the number of surface points of the IMN located outside MC to optimize the registration. If a point on the IMN falls outside the MC, it means that the IMN squeezes the inner wall of the femur there. Here are the processing steps:

(1) Acquisition of the surface points of the IMN falling between the adjacent control points of the MC centerline: As shown in Figure 3, A and B denote two adjacent control points on the MC centerline. Plane D ₁ and plane D ₂ are cross sections that pass through point A and point B respectively and are perpendicular to the line connecting points A and B. All the points between planes D ₁ and D ₂ in the transformed IMN point cloud are saved in point set P _middle .

We find that some points fall outside the MC after calculation because the isthmus length is too short to take into account the registration of both ends. The best registration position can be found by changing the length of the isthmus centerline and registering many times to reduce the number of outside points. However, there are many options for the growth direction and length of the isthmus centerline. If all the options are performed separately, the calculation efficiency will be greatly reduced. Therefore, we propose an iterative adaptive registration method to quickly find the effective length and position of the isthmus centerline. In this method, the isthmus centerline can determine the growth direction and growth speed autonomously according to the distribution of outside points. The MC is divided into upper and lower parts bounded by the midpoint of the isthmus centerline. We assume that the number of totally outside points of the IMN is N o u t p , and the number of points beyond the upper MC is N o u t u p . If N o u t p ≥ N o u t u p > 4/5 N o u t p , the centerline of the MC isthmus increases 1 mm upward at a time (Figure 4A). If 4/5 N o u t p ≥ N o u t u p > 3/5 N o u t p , the isthmus centerline increases 2 mm upward and 1 mm downward at a time (Figure 4B). If 3/5 N o u t p ≥ N o u t u p > 2/5 N o u t p , the isthmus centerline increases 1 mm upward and downward respectively each time (Figure 4C). The isthmus centerline increases 1 mm upward and 2 mm downward if 2/5 N o u t p ≥ N o u t u p > 1/5 N o u t p (Figure 4D), and 1 mm downward only every time if 1/5 N o u t p ≥ N o u t u p ≥ 0 (Figure 4E). The growing will stop when the MC isthmus centerline grows to the endpoint. The growth direction and growth speed are determined according to the proportion of N o u t u p to N o u t p . The proportion is obtained by a large number of experiments. We adopt the corresponding growth scheme when N o u t u p is in a certain range, and the IMN and MC can quickly achieve the best registration effect.

When the isthmus centerline increases, the centerline with the same length at the corresponding position of IMN also increases. According to the determined growth scheme, each increase of the isthmus centerline is re-registered with the IMN centerline to obtain a new transformation matrix T. The IMN is registered with the MC based on T. The number of outside points N o u t p , the mean exceeded distance D o u t p , and the mean distance (marked as E c o n t r o l p ) under the registration are calculated and recorded. The value of E c o n t r o l p is calculated as Eq. 5: E c o n t r o l p = ∑ i = 1 n d i / n (5) where n is the number of control points in the IMN centerline. d i is the Euclidean distance between ith control point of the MC centerline and ith control point of the IMN centerline. The normalization N o u t M of N o u t p is processed as follows: N o u t M j = N o u t p j − N o u t p Min N o u t p Max − N o u t p Min (6) where N o u t p j is the number of outside points under the registration after jth increase. N o u t p Max and N o u t p Min denote the maximum value and minimum value of all values of N o u t p after each increase. The same operation is performed on D o u t p and E c o n t r o l p to obtain D o u t M and E c o n t r o l M . The fitness of the IMN registered with the MC after jth increase is calculated as Eq. 7: Ψ j = o u t M N j 2 + o u t M D j 4 + c o n t r o l M E j 4 (7)

The optimal growth times can be found by comparing the values of Ψ of each increase. At the same time, we notice that the selection of the initial length of the isthmus centerline will also affect the accuracy of registration. The initial length is set to 70 mm, 80 mm, 100 mm, and 110 mm respectively. The isthmus centerline increases again according to the determined growth scheme and registers with the IMN centerline multiple times. Finally, the optimal growth times and transformation matrix under the optimal isthmus length are obtained. The position of IMN based on the transformation matrix is the optimal registration position (Figure 1F). The main computing process is shown in Figure 5.

The correct location of the entry point is very important for smooth insertion of the IMN and no deformation after insertion. In this study, a new solution is developed to choose an optimal entry point. After the IMN is registered at the optimal position, ten continuous control points are selected near the upper end of the IMN centerline (Figure 2F). The ten control points are fitted as a straight line, and the line is extended to the surface of femoral greater trochanter. Distances between all the proximal femoral surface points and the line are calculated successively. The point with the shortest distance is adopted as the optimal entry point. In clinical practice, the femoral image coordinate system and the patient coordinate system will be registered. Surgical instruments equipped with spatial trackers can locate the actual entry point on the femur according to the determined entry point on the femoral image. Then, the opening operation of the entry point can be guided by spatial navigation systems.

In clinical surgeries, the first thing is to judge whether the length of an IMN is suitable or not. Ideally, the distal end of the IMN should be located between the patella upper edge and the bone scale line. Too long IMN will destroy the cortical bone of intercondylar fossa or hinder incision suture. Too short will affect the stability of fracture fixation and the function of hip joint and knee joint, and may even lead to postoperative re-fracture. In this study, an appropriate length of IMNs is determined between the length of MC centerline −3.5 cm and the length of MC centerline. Similarly, the diameters of IMNs should also be considered. A large diameter will lead to implantation failure, and a small one will lead to insufficient support strength. In this study, the maximum diameter of the lower half of IMN is determined between the minimum cross-section diameter of MC -3 mm and the minimum cross-section diameter of MC. If an IMN does not meet the length and diameter requirements, no further processing will be carried out. IMNs that meet the length and diameter requirements perform steps 1 to 5 of the method in sequence, and the optimal registration positions of all the IMNs are obtained. The number of outside points N o u t p , the mean exceeded distance D o u t p and the mean distance E c o n t r o l p of each IMN under their corresponding optimal registration are calculated and recorded. The support strength of IMN (marked as S I M N d ) is also an influencing factor that must be considered, which involves the early weight-bearing exercise of patients. The support strength of IMN is mainly related to the anti-torsion strength and the anti-bending strength, which are proportional to the cube of IMN diameter. We use the cube of the difference between the minimum cross-sectional diameter of MC and the maximum diameter of the lower half of IMN to characterize the support strength, and the maximum difference is set to 3 mm. The universal measure N o u t U M of N o u t p is processed as follows: N o u t U M j = N o u t p j − N o u t p min N o u t p max − N o u t p min (8) where N o u t p j is the number of outside points of jth IMN under its optimal registration with the MC. The minimum of N o u t p donated N o u t p min is 0. N o u t p max represents the maximum number of IMN surface points that are allowed to exceed outside the MC. According to the registration results of multiple IMNs and medullary cavities, N o u t p max is 120 which is an empirical value obtained from many tests. Similarly, the registration results also get that the maximum value of D o u t p is 0.47 mm, and the maximum value of E c o n t r o l p is 5 mm. Because of the relatively small sample size, all the maximum values are temptative values which require further studies. All the minimum values of N o u t p , D o u t p , E c o n t r o l p and S I M N d are 0. The same operation in calculation of N o u t U M is performed on D o u t p , E c o n t r o l p and S I M N d to obtain D o u t U M , E c o n t r o l U M and S I M N U M . The suitability Φ of jth IMN for the MC is calculated as Eq. 9: Φ j = 2 × N o u t U M j 5 + D o u t U M j 5 + E c o n t r o l U M j 5 + S I M N U M j 5 (9)

The smaller the Φ is, the higher matching between the IMN and the MC is. The most suitable IMN for the patient is determined by comparing the values of Φ (Figure 1G).

The effectiveness of the registration between IMNs and MCs based on the centerline growth scheme of MC isthmus was verified in this experiment. Patient 1 was taken as an example. The No. Two IMN was the most suitable IMN for patient 1 determined by the proposed method. The adopted growth scheme of the isthmus centerline was to increase 1 mm upward and 0 downward each time. According to the calculation result, the best registration length was reached in 27 growth times. Five different growth times were selected according to the rule of sequential increase, and there were two growth times before and after the optimal growth times respectively. The five growth times of patient 1 were first, 14th, 27th, 40th and 53rd growth. The five growth times were denoted as times 1, 2, 3, 4 and 5. Similarly, the five growth times of patient 2 were first, 13th, 27th, 39th and 51st growth. The five growth times of patient 3 were first, 13th, 25th, 37th and 49th growth. While patient 4 were first, fifth, eighth, 20th and 35th growth. The values of N o u t p , D o u t p and E c o n t r o l p of registration between the MC of patient 1 and its optimal IMN at 5 different growth times were calculated and recorded in Table 1. The same operation was performed on the remaining patients respectively, and the values of N o u t p , D o u t p , and E c o n t r o l p were also recorded in Table 1.

As shown in Table 1, all N o u t p , D o u t p and E c o n t r o l p decrease at first and then increase with the increase of growth times, and all get the minimum value at the optimal growth times. Therefore, the registration is affected by the position and length of the isthmus centerline. The result indicates that the proposed method can obtain an appropriate centerline length of the MC isthmus.

The proposed method was verified furtherly by comparing with a mature software (Geomagic Control 2014). According to the calculation of the proposed method, the most suitable IMNs of the four patients are the No.2, No.3, No.1 and No.3 IMNs successively. Patient 1 was taken as an example. The No.1, 2, 3 and 4 IMNs met the length and diameter requirements of patient 1. Among them, the No.2 IMN is determined by the proposed method, and the other nails were selected as a control group named contrast 1, 2 and 3. The femur model of patient 1 and the No.2 IMN model after transformation were imported into the software. The Boolean operation was employed to obtain an intersection from the surface of the MC model to the surface of the IMN model. The volume and wall thickness of the intersecting were calculated and recorded. The other three IMNs under their corresponding best transformation and the femur of patient 1 were imported into the software respectively. The volume and wall thickness of the intersecting parts were calculated. The remaining 3 patients performed the same operations. The numbers of the outside points of IMNs calculated by the proposed method were compared with the intersection volume in Figure 7A. The wall thickness distribution of the intersecting parts was shown in Figure 8. Next, the best registration position was verified. Patient 1 and its most suitable IMN were taken as an example. The best registration position obtained by the proposed method was recorded as position 1. The other four positions are selected: An optimal position calculated by the isthmus initial length which is different from the optimal growth scheme was arbitrarily selected as position 2. An optimal position different from the growth direction and growth speed of the optimal growth scheme was arbitrarily selected as position 3, and any two positions different from the best position under the optimal growth scheme were selected as position 4 and position 5. After being transferred to the five positions in turn, the IMN model and the femur model were imported into the software. The volume and wall thickness of the intersecting parts were calculated and recorded respectively. The remaining 3 patients performed the same operations. The numbers of the outside points of IMNs in five registration positions were compared with the intersection volume calculated by the analysis software in Figure 7B. The wall thickness distribution of the intersecting parts was shown in Figure 9.

As shown in the result, the volume and wall thickness of intersecting part of the IMN determined by the proposed method are smaller than those of the other three IMNs in the control group. The number of outside points calculated by the proposed method is proportional to the intersection volume. The same result is obtained for the verification of the best registration position. The volume and wall thickness of intersecting part in the best registration position are smaller than those in the other four positions. For some IMNs, it is found that the number of outside points is zero while the intersection volume is not zero. This is because the MC is not a regular cylindrical tube, and a slight error occurs in the operation of the layered synthesis of circles. However, the wall thickness caused by the error is still small. The average intersection wall thickness of the four patients under the best registration is 0.191 mm. This accuracy satisfies the clinical requirements.