Designing of Patient-Specific Implant by Using Subdivision Surface Shaped on Parametrized Cloud of Points

Patient-specific orthopaedic implant application improves the recovery of the fractured bone in comparison to the conventional implant of standardized shape and size. The main challenge regarding creation of this kind of implant is its design process because the conventional modelling techniques based on using NURBS for shaping the unique bio-forms take a lot of time and effort. The research found out that the application of combination of two design techniques may accelerate the whole design process for more than 60% comparing to the conventional approach while improving the geometric congruency of the implant and the bone. The first is about digital reconstruction of the bone geometry by the polygonization of the synthetic cloud of points built on radiographic images of the injured bone. The second is applying semiautomated surface subdivision modelling technique for shaping the implant. The research was done for the case of designing the internal fixator of Mitkovic type aimed for lateral tibial plateau fracture.


INTRODUCTION
There are three types of orthopaedic implants: endoprostheses, scaffolds and fixating assemblies (which are usually consisted of plates, screws, pins etc.). Personalizing or customizing a bone implant to the specific patient comprises adapting implant's design to the patient's anatomy (i.e. geometry), harmonizing implant's mechanical properties to the patient's bone mechanical properties, and finally, making the implant biocompatible to the patient-specific biochemical complex. Besides these three major aspects of a bone implant personalization, there are additional two constraints that certainly should be taken into consideration whilst designing the implant specifically to a particular patient: 1) Implantation procedure combined with a surgical operation conditions that may be specific for each patient, 2) Manufacturing technology that can be applied for the implant fabrication.
In recent times, emerging manufacturing technologies (especially additive manufacturing technologies -AMTs) make the application of non-standard, patient-specific implants (PSIs) to become more affordable. AMTs enable fabrication of very complex shapes in a short time which makes them ideal for PSI manufacturing. However, in the case of PSI production process, due to their complex and unique geometry, the designing of PSI turns up as the bottleneck phase, which takes most of the time and effort. Since it seems there is not a precisely defined efficient design procedure that is proposed for the PSI modelling except for some specific cases [1][2][3][4][5] the challenge was to identify the most efficient available CAD technique that may be applied.

Patient-specific Implant Design Process and Related Research
Besides a few patents which describe the design method of patient-specific orthopedic implants with more details [1,2], there are not many published researches 1 In context of orthopaedic surgery, region of the specific interest -RSI denotes a specific anatomical section defined as part of the bone where (papers) that were focused on PSI modeling. In the patent [1] the core of invention is described within two methods for creating accurate CAD models of the interface surface of the implant for an example of creating distal femoral implants for partial or total knee replacement are presented. Both methods are related to the tensor B-Spline Surface fitting using least-squares-based techniques. These methods are iterative whose goal is to determine the control net size as well as the desired smoothing parameter so that the maximum error bound is satisfied. The difference between these two methods is whether the implant interface surface is being constructed as one-or multi-NURBS patch. However, in these patents, which is after all to be expected, the details of the research are not presented, the choices for the presented method are not argued in detail, nor the challenge of refining, smoothing and connecting several NURBS patches into a whole was analysed.
The design process of patient-specific implant is usually conducted in two phases, as it is shown in diagrams ( Fig. 1, Fig. 2): firstly, by digital reconstruction of the patient's bone geometry and, later, by the modelling of corresponding implant (fixating plate).
In the first phase, the objective is to build the tessellated surface model of the bone or just an area of ita region of the specific interest (RSI) 1 . Usually, in real cases, the radiographic images (X-ray, CT scans, or less often MRI) show the fractured bone, which by fortune may be without dislocation and with all pieces of the bone in place. In that case, the cloud of points (CoP) made from the radiographic images can be used directly. In the case fractured bone pieces are dislocated or even crushed (so they are lacking), it is very difficult to reconstruct the bone shape before it has been injured. For this case, it is possible to apply two approaches: the first is to use images of the healthy bone that is mate to the one that is injured (for example, if the left tibia is fractured, to use the right one, which is not fractured). CoP of the uninjured bone are the implant should be placed. At the same time, from the CAD perspective, RSI refers to a part of the cloud of points that should be used for surface formation. being mirror-transformed so the geometry of the injured bone is a mirror of the healthy one. The second approach is to generate so-called synthetic cloud of points (S.CoP). The part of the CoP (which is missing) is synthetized according to the part of the CoP which remained undamaged. The construction of S.CoP is most often performed on the basis of statistical analysis of mutual geometric relations of points in the cloud obtained from numerous CT scans of healthy (same) bones [6][7][8][9][10]. The S.CoP generated in this way are often of reduced geometric accuracy, but are cleared of the errors common to CoP obtained directly by CT scans (real CoP). Whatever CoP (synthetic or real) is used, the following three design operations are required to make a usable tessellated model of the bone surface: 1) refining the CoP, tessellationcreating the crudely tessellated surface model of the bone over CoP, 3) healing and smoothing (editing) the crudely tessellated surface.
Although the digitally reconstructed bone surface in the format of a tessellated (mosaic/polygonal) model can be very accurate in geometric terms (as in the case where real CoP can be used), it is unsuitable for further manipulation of the geometry. In the paper [11], а method of processing the CT scans (radiographic images) for restoration of cranial defects is presented. The open source toolkit VTK is used for upgrading and applying existing algorithms for radiographic medical images processing. The outer surface of the reconstructed model through Delaunay triangulation while the contour clipping is accomplished by using a VTK class vtkClipPolydata. This kind of approach is shown as very useful and efficient if the tessellated surface can be directly used for the implant fabrication by some AMT like in the case for a specific cranial endoprosthesis presented in the paper. However, if it is necessary to design a completely new implant that should partly resemble the tessellated surface of the bone, as in the case of the fixation plate, then it is necessary to significantly intervene and change the geometry. In this case, it is necessary to create manageable surface model of the bone or RSI.
Another group of papers whose theme is close to this one, though not related to the design of PSI, addresses the problem of converting 2D radiographic (X-ray) images to 3D models of bone [12][13][14]. The last one [14] is based on Laplacian surface deformation method and it may be very useful for the design of implant interface surface, especially in the case when there are just 2D x-ray images on disposal. However, these procedures also result in only tessellated surface of the bone, i.e. a 3D surface assembled of numerous triangles or polygons that are difficult to change further.
In the second phase, the creation of manageable wrapping surface model of the patient's bone or RSI is necessary because this model serves as a referential geometric model for designing the interface surface of the implant, and finally for the design of implant itself. Creating the manageable surface model of the bone can be performed through three different methods (Fig. 2). Associating the implant's with bone's geometry The first method (option 1 in Fig. 2) relies on applying conventional NURBS modelling techniques: forming the wrapping NURBS surface by dragging over multiple crosssection curves of the tessellated model. The second method (option 2 in Fig. 2) employs automatic forming and merging of numerous small and separate iso-parametric NURBS patches shaped per the tessellated model into the one wrapping NURBS surface. There are a series of papers of a group of authors [3, 4, 8, 9, 15-18, 20, 21] which describe these two methods in detail. It may be said that these methods are the current standard for design of patient-specific implant. The main shortcoming of the first method is about considerable additional effort needed to refine, smooth out and finally merge separate NURBS patches. Despite the considerable improvement comparing to the first method in terms of time and effort needed to make it, the model of manageable surface built by the second method is difficult to modify in further steps, which is usually important for the additional implant design adaptation. Unifying the patches into one smooth and connected NURBS can take additional time and effort (for example, for the case of tibia fixating plate, merging and smoothing the mosaic NURBS took an average of 72 minutes).
The third option (option 3 in Fig. 2) is to apply the technique of sub-division and shaping Non-Uniform Recursive Subdivision Surface (NURSSs) primitive, or even better Non-Uniform Rational Catmull-Clark Surfaces with T-junctions (T-NURCCs), to match the tessellated model. Shaping is mostly performed by collocating the vertices and sides of the control polyhedron of a built SubD surface. At the same time, by creating the SubD surface that mimics the specific bone region, the congruent geometry of the laying bone-implant interface surface (BIIS) of the corresponding implant is being designed. This way-built surface is ready for further manipulation, that is, for designing the corresponding implant. The model made as SubD surface may be changed and adjusted easier than the model made by first two mentioned methods, but the geometric accuracy of small, local details such as crisps and small fossas is hard to achieve. However, in the context of orthopaedic surgery and implant congruency with the bone surface, the geometric accuracy of the bone model made as SubD surface is quite sufficient. The main challenge for this kind of approach is related to the automation of sculpting the NURSS or T-NURCCS surface by collocating the control polyhedron points. The newest and very promising solution is related to the usage of T-NURCC and NURSS [5,15,16] for shaping the onepatch surface in accordance with the tessellated surface of the bone. This paper also belongs to the group of papers which recommend subdivision surface techniques on NURRS and T-NURCC surfaces.

PROBLEM STATEMENTAND THESIS
The goal of the research is to determine the most efficient (available) CAD technique that should be applied for the case of designing the PSI. Analysing all the operations in the patient-specific implant design process, two were identified as the ones whose improvement may bring out the greatest acceleration to the whole process.
The first is related to the generation of the cloud of points based on radiographic images of the injured bone, and the thesis is that one should use the synthetic (parametric) cloud of points. Generating the synthetic (parametric) cloud of points formed by parametric functions, whose arguments are morphometric magnitudes related to the prominent referential anatomical entities, enables fast creation of tessellated surface model of the bone (as if it has not been injured) that is simplified, but still quite sufficiently congruent with the bone geometry and it is error-free.
The second is related to the operation of designing the bone-implant interface surface as a basis for further patient-specific implant design, and the thesis is that one should apply surface subdivision modelling technique instead of conventional techniques of NURBS creation for shaping the unique bio-forms. Surface subdivision and sculpting enables fast and accurate shaping of a NURSS (T-NURCCS) primitive to the required bio-form, avoiding generation of cracks, gaps, tweaks, and wrinkles in contrast to the techniques that are based on NURBS patches modelling.
The paper presents the design method that may be considered as a very efficient regarding the PSI modelling. It comprises application of in-house developed method for radiographic images (X-ray, CT and MRI) processing [10], upgraded to harness the modelling technique based on surface sub-division. The application case was worked out for designing the internal fixator for lateral tibial plateau of Mitkovic type [21]. In addition to the approach presented in [21], the CAD technique of surface sub-division (SubD) is improved by introducing the control based on usage of prominent morphologic points of the bone as shapedragging vertices. In addition, it was already explored that SubD gives the advantages over conventional NURBS modelling in the context of bio-shapes design [15,16].

DESIGN CASE: INTERNAL FIXATOR OF MITKOVIC TYPE AIMED AT LATERAL TIBIAL PLATEAU FRACTURE
The research was conducted for the case of designing the fixating plate aimed for the recovery of human tibia injuries in the proximal areas, close to tibial plateau. The typical fracture in the proximal tibia zone that indicates the application of internal fixator for lateral tibial plateau is showed in Fig. 4a. For the Mitkovic type of this class of fixating plate the proper position is defined in [22] and presented in Fig. 4b. Tibial plateau fractures are complex injuries produced by high-or low-energy trauma. They imply cracks or breaks in the upper portion of the tibia which may or may not involve the knee joint. They are also known as proximal tibial fractures. Tibial plateau fractures are often associated with serious soft tissue damage. They account for 5-8% of the lower limb fractures and represent about 1% of all fractures [23,24]. Bone plates are commonly used in internal fixation of bone fractures. Compression plates are type of plates which are designed to apply compression to fractures to provide stability. They are made in different sizes and shapes so they can be used for different patients. Internal bone fracture fixation includes application of compression plates [25]. This group of implants can be separated on Dynamic Compression Plates (DCP) [26] and Locking Compression Plates (LCP) [27], or some type of their combination. With DCP direct contact between the bone outer tissue (periosteum) and the implant is required. The pressure on the periosteum can lead to devitalisation of the underlying bone and it can result in plate breaking [25]. For LCP, screws are inserted into bone and locked into the plate, so that they provide more stability to the assembly of the bone and the implant. There is no requirement for any direct contact of the bone surface and plate implant, but it is better if the implant geometry and shape follow the bone geometry [25][26][27][28][29]. According to the references related to the biomechanical properties of unique plates (LCP) vs. standard plates (DCP) [22,29], it is concluded that the unique plates provide better fixation and they can withstand more load.

Creation of Synthetic Cloud of Points
The procedure for creating the synthetic cloud of points (S.CoP.csv) described in [30] and showed in Fig. 1 has been used to prepare the tessellated surface model of the human bone, that is, in this particular case, for the model of proximal tibia. For the determination of relations between the referential geometric entities (like morphologic prominent points and axes) of tibia and all other points of CoP that are used for the tessellation, the collection of radiographic images of healthy tibia samples of University Clinic of Nis was used. The morphometric measurements that were carried out over each tibial CT image sample provided extensive datasets representing the spatial relationships between the referential morphometric parameters (RMPs) and other points of the bone enveloping surfaces. Performing the statistical analysis (multilinear regression [10,31], and in some cases neural networks [32]) on these datasets, the parametric equations, which define multivariable linear relations between RMPs and coordinates of all other anatomical points, are derived, Eqs.
(1) to (3) These parametric equations could be comprehended as a sort of generic S.CoP for a human tibia. To generate a S.CoP for tibia of a specific patient, it is enough to capture the magnitudes of RMPs from the radiographic images of the patient's tibia and insert them into the corresponding parametric equations. A Python shell application, which runs the set of parametric equations, enables designer or surgeon to enter the magnitudes of RMPs (in format of *.csv file) and get the coordinates of all anatomical points of the bone or RSI of the bone (also written in format of *.csv file).The greatest advantage of using the procedure for getting the S.CoP is that majority (a sufficient set) of RMPs, if not all, may be taken from the CT scan of injured, fractured bone. Considering that S.CoP actually represents how the injured bone approximately should look like without (i.e. before) the fracture, it will serve as a basis for tessellation of the boundary surface model of the bone. The tessellated surface model of the bone or RSI of the bone is free of errors, that is, gaps or non-manifold vertices. Even though, this model is not too accurate in geometric terms, it is suitable for building the NURBS or NURSS (T-NURCCS) over it.

Creation of a Manageable Surface Model
In a case the NURBS patches are going to be used, the procedure requires tibia tessellated model to separate into patches by using anatomical points (landmarks) which belong to the different regions (areas) of the bone geometry (Fig. 3). The NURBS boundaries pass through these boundary anatomical points (landmarks). These regions are usually defined arbitrary, because division into sections is being done in accordance with the surgeon's recommendations following the anatomical terms and definitions as presented in [22,33] and in Fig. 4. The anatomical regions within the tessellated model of the bone are represented by a set of points from S.CoP as presented bellow in Eqs. (4) to (7). ARk : k th anatomical region; Pk,i : i th Point in the k th anatomical region; V : Set of anatomical regions

 
In the next step of the procedure the enveloping surface of each anatomical region is being formed by running the Python script. The script utilizes the corresponding Python libraries and points' coordinates to generate the boundary curve and NURBS surface within it.
However, the regions of specific interest (in context of orthopedic surgery operation for the specific patient) usually match the anatomical regions partially. Moreover, if we are going to design patient-specific fixation plate it is easy to notice that the projection of the fixation plate crosses several multiple trimmed NURBS patches that build enveloping surface of the bone geometry (Fig. 4b,  Fig. 5 and Fig. 6). If there are (and very often this is the case) some cracks or small gaps between the NURBS patches, the integrity of the fixation plate projection curve onto NURBS patches, can be corrupted [15,35]. Therefore, the interface surface of the implant, the boundary of which should be defined by the fixation plate projection curve onto NURBS patches, could not be joined. Even though this can be overcome by precise modelling and using mutual B-spline as demonstrated in [15], that is usually a complex task and it requires a lot of time, which is in a contrast to the required quick response in clinical cases.
NURSS, owing to its ability to follow the complex topology (often present in bio-shapes) and its easiness to use, has become an attractive alternative to NURBS. The subdivision representations have a great advantage that can be seen for the case where two subdivision surfaces share a boundary edge in base meshes. In that case, both SubD surfaces contain exactly that piece of boundary curve. This makes it possible to apply the same trimming curves on two intersecting surfaces, and thus, providing the gap-free models [15], with a curvature continuity and topology accuracy preserved.

Bone-implant Interface Surface Designing
As it was stated herein before, for the case of designing the patient-specific fixation plate (in this particular case the plate of Mitkovic type), the bone-implant interface surface should be modelled firstly. Mitkovic type of internal fixation for lateral tibial plateau is consisted of three distinctive sections (Fig. 5): the top section that should rest on the proximal lateral part of the tibia, the lower section that should rest on the mid-or distal part of the tibia shaft, and the mid-section that connects these two fulcrum sections.
In the projection showed in Fig. 5, the top section of the plate is characterized by spatially curved flank that should follow the anatomy of patient's bone in its RSI (area of proximal-lateral zone just below the tibial plateau). The fulcrum sections (top and lower section) are additionally characterized by the bumps on the fixator side towards the bone, which round the holes for fixating screws placement. These bumps ensure the required clearance between the plate and the bone that is crucial for unclogged blood flow and nutrition to the periosteum [10,[36][37][38] of the RSI. The silhouette curve (initial version) of a top section of the Mitkovic fixating plate projected on the tessellated surface of proximal lateral area of tibia altogether with the defined anatomical points that are going to be used for building the SubD interface surface are showed in Fig. 6.  The first step of creating the SubD surface is to model the silhouette (contour) of the fixating plate in a plane which will define the direction of the projection onto the tessellated bone surface. This plane is being placed arbitrary. However, it should be parallel "as much as it is possible" with the flattish lateral tibial shaft area (Fig. 7).
The projected contour of the fixating plate is created onto the tessellated bone surface, as presented in Fig. 8. The initial dimensions are measured from real plate implant of Mitkovic type (Fig. 8). The base uniform quad mesh for the desired subdivision surface is being constructed on the tangent (projection) plane and defined boundary (Fig. 9). The control points' positions of the base quad mesh are mapped according to the pre-defined limits (or number) of subdivision sections. The created subdivision area, i.e., "polygonal net", intersects with plate boundary curve and forms the limiting SubD surface mesh, which will be projected onto the tessellated model of the bone, by following the projecting direction (Fig. 8), in the following step. It is important to note that overall width and length dimensions of the quad mesh are being set a bit larger than the original plate size to enable additional mesh adaptation if required, (Fig. 9). The next step (as stated) is projecting the initial subdivision surface mesh to the tessellated surface as presented in Fig. 10. The SubD surface is constrained and controlled by the vertices of the polygon that have got associated to the tessellated surface built over synthetic CoP. The integrity of the mesh has been checked in Blender. If there is a need to intervene in order to form a valid, one-patch mesh, it is possible to make corresponding corrections and then transfer it back to CAD software (CATIA V5 R25, Fusion360).
In order to automate the creation of SubD surfaces, Python library -PyMesh [15] can be used, which provides all tools necessary for creating valid subdivision surface mesh, whose vertices are disposed on the sketched boundary curve of the fixation plate in a predefined manner. In addition, if there is a need, this library can be enriched by the very useful recently presented algorithms that can help accomplishing these sequences [5,39].
In this way designed subdivision surface is ready for further manipulation in terms of overall and local shaping, as well as for mesh refinement (Fig. 11). Currently, this operation is being performed in CATIA (Blender is also an option), but in the future work, Python script will be developed to automate this step. The NURSS offset from bone-implant interface surface is permanently associated with it so each change on BIIS is being carried pass onto design of the NURSS surface which is at the same time the base for the implant model (Fig. 12).

Figure12
Offsetting the interface surface of the implant for the required magnitude (e.g. bumps' height) creates the PSI base The last series of operations are related to creation and local shape adjustment of initial solid model of patientspecific implant (Fig. 13). However, all further design changes and refinement of interface surface of implant that are going to be made, will be entangled with PSI. Figure 13 Thickening the offset NURSS surface forms the initial shape of PSI

CAD Techniques Comparison
The result of the applied approach can be summarized in the following observations: 1. Two working hours are required for designing the initial shape of PSI of such a kind from the moment of importing S.CoP; a. With automation of shaping the interface surface of implant to the tessellated surface of the bone by controlling the vertices of the NURSS polyhedron through a Phyton script, the whole process may be shortened even more, improving the matching at the same time; 2. The model is fully associated with tessellated surface model of the bone; 3. Customization of the implant shape overall or locally is quite simple, but it requires some talent for digital sculpting. Compared to two NURBS modelling techniques, this approach shows obvious advantages in terms of productivity, quality, and flexibility (customization), (Fig.  14, Fig. 15, Fig. 16, Fig. 17).
The productivity for the different CAD techniques is measured by the time needed to model interface surface of PSI with the same CAD software (CATIA, Fusion 360), engaging the same CAD skills and hardware equipment. The test has been repeated for eight tibia samples -cases of similar fracture of tibial plateau (Fig. 14, Fig. 16, Fig. 17).
The quality of these CAD methods has been measured by two magnitudes: 1. By the time needed to achieve target accuracy (congruency) which is defined as deviation between the bone tessellated surface built on the real cloud of points and bone-implant interface surface that is equal or lower than 0.5 mm 2. By the time needed to refine and smooth the model removing the errors like cracks, gaps, and wrinkles (if any), (Fig. 14, Fig. 16). The flexibility is measured by the time needed to perform additional design changes of the bone-implant interface surface after its initial creation (and consequently all the associated changes of the implant design). For the test purpose, additional shaping and widening of the boneimplant interface surface ends have been selected as the reference (Fig. 14, Fig. 16).

Figure 16
Comparison of averaged times needed to accomplish three referential tasks by three different CAD techniques. For the sake of comparison, the design operation time is normalized regarding the slowest method.

Figure 17
In case the CAD technique which creates mosaic NURBS patches is taken as a reference with a productivity of 1, the productivity of the other two methods can be expressed relatively.

Limitations and shortcomings
There are two main shortcomings of the proposed method (procedure): The first is related to creation of synthetic cloud of points -it is mathematical, and not the real one. It is built on learned spatial relations between referential morphometric parameters and the prominent points in the clouds of series of similar healthy bone shapes. The accuracy of CoP built in this way may vary a lot. That can affect the accuracy of the implant geometry substantially. In future research, the solution could be explored in building the computer-aided program module for automated re-assembling the real cloud of points generated from the radiographic images of the injured (dislocated and fractured) bone. In this way, one may expect to get the cloud of points that represents the real bone, but before it has been injured. However, it should not be forgotten that the real cloud of points abounds in many erroneous details that make the tessellated model of the bone surface corrupted and whose healing requires great effort and time.
The second is related to usage of techniques for surface subdivision and morphing up to a required bioform. For the efficient shaping of this kind of surfaces it is expected from the CAD expert to be able to conceive complex spatial relations and shapes, i.e. a bit of talent for digital sculpting is needed, which is a more than usual request for a mechanical designer to fulfill.

CONCLUSION
The request for designing the personalized internal fixator of Mitkovic type aimed for lateral tibial plateau fracture, for the specific patient, has induced search for the most efficient design method that could be applied. As a result, the research has advanced the method of anatomical features (MAF) that generates so-called synthetic cloud of points and, consequently, the tessellated surface model of the bone that is sufficiently congruent with the real geometry of the injured bone. Also, the research gave birth to a procedure for creation and trimming of quad mesh of a NURSS (open) planar primitive bounded by the contour of the planar projection of the patient-specific fixator. Finally, the conducted research pointed out a possible method for the automatic shaping of a NURSS primitive over the tessellated surface model of the bone by mapping the control points of the NURSS polyhedron. With appropriate additional local subdivisions and shaping interventions of the SubD surface entities such as translation, rotation, scaling, affine transformation (points, edges, and faces), it is possible to highly accurately fit the implant surface interface to the bone surface in a fast and easy way. The bone-implant interface surface made in this way becomes a base for creation of patient-specific (personalized) internal fixating plates of any kind, either dynamic compression plate (DCP) and locking compression plates (LCP), with required geometrical accuracy and anatomical correctness.
The most important, with this method orthopaedic surgeon has got the tool to efficiently adapt the implant geometry in accordance with the patient anatomy and mechanical properties required for the patient's specific mechanical load case that is designed for the recovery process. Besides its basic purpose -production and implantation of personalized plate implants and fixators, 3D models of implants created in this way can be used for pre-operative planning, surgical guidance, and educational purposes in medical science.