Surface Reconstruction from Parallel Curves with Application to Parietal Bone Fracture Reconstruction

Maxillofacial trauma are common, secondary to road traffic accident, sports injury, falls and require sophisticated radiological imaging to precisely diagnose. A direct surgical reconstruction is complex and require clinical expertise. Bio-modelling helps in reconstructing surface model from 2D contours. In this manuscript we have constructed the 3D surface using 2D Computerized Tomography (CT) scan contours. The fracture part of the cranial vault are reconstructed using GC1 rational cubic Ball curve with three free parameters, later the 2D contours are flipped into 3D with equidistant z component. The constructed surface is represented by contours blending interpolant. At the end of this manuscript a case report of parietal bone fracture is also illustrated by employing this method with a Graphical User Interface (GUI) illustration.


Introduction
Craniofacial region is a complex anatomical part made up of various bones joined together. Fig 1 illustrates the various bones that make up the human skull. Craniofacial fractures occurs due to various etiological factors like road traffic accident (RTA), sports injury falls etc. Various diagnostic tools like X rays, Computerized Tomography (CT) scans, Magnetic Resonance Image (MRI) have been used to diagnose the craniofacial fractures. Since, craniofacial fractures do not follow a specific pattern and lately, emerging virtual reconstruction technologies opened new avenues for mathematicians, physicists and software engineers to reconstruct the fracture defects. Already established approaches for implant design is Computer Aided Design (CAD)/Computer Aided Manufacturing (CAM) process chain [1]. Other alternative methods that design implant without CAD process are mirroring [2] and surface interpolation also called as deformation [3].
The methods used for the construction of fracture segment are either anatomical or mathematical features. The mirroring method is useful for the fracture on one side of the skull only and use anatomical features. While mathematical constraints are used in interpolation and deformation methods to construct the fractured part from the non fractured part of skull. Reference model is also used to construct the fractured part of skull [4]. For more related work (see [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] and references therein).
Contours blending function have been used to construct the 3D parietal bone fracture. Since Digital Imaging and Communications in Medicine (DICOM) data used in our work are in 2D form, so, before constructing the 3D surface, firstly, we have constructed fracture part curves using geometric continuity (GC 1 ) rational cubic Ball curves with three free parameters. Genetic Algorithm (GA) is used to optimize the free parameters. Then, we converted the 2D contour curves in 3D form by taking equidistant z component and then construct the 3D parietal bone fracture.
We present a case report of parietal bone fracture to show the applicability of proposed algorithm. The Dicom data are obtained from Hospital Universiti Sains Malaysia (HUSM). We employed matlab for both the programming and to develop Graphical User Interface (GUI) for parietal bone fracture reconstruction using proposed algorithm. Surgeons can use GUI for fracture reconstruction without having in depth knowledge of its mathematical aspect.
The introduction is followed by the representation of Ball basis functions and rational cubic Ball curve. It is followed by an explanation on contour blending, boundary extraction, parametrization, Normalized mean squares error. After this, we have explained the Graphical User Interface (GUI). An algorithm on 3D parietal bone fracture reconstruction using rational Ball curve and contour blending is proposed after the explanation of GUI.
• Partition of Unity: Ball basis functions form a partition of unity, it means that sum of Ball basis functions will be 1.
S i ðyÞ ¼ 1: Theorem 2: Let P i be the set of control points and S i (θ), i = 0, . . ., 3 are Ball basis functions defined in Eq (1). The Ball curve sðyÞ ¼ P 3 i¼0 P i S i ðyÞ have the following properties • Coordinate system independence: As the Ball basis functions form partition of unity, so the Ball basis curve will be coordinate system independence. It means that by changing the coordinate system of control points curve will remain same.
• Convex Hull Property: As the Ball curve obeys the coordinate system independence and Ball basis functions are all non negative. So Ball curves will obey the convex hull property. It means that curves formed by Ball basis always lie within the convex hull of their control points. • Endpoint Interpolation: The Ball curves always passes through the first and last control points.
Theorem 3: Let P i be the set of control points and S i (θ), i = 0, . . ., 3 are Ball basis functions defined in Eq (1). The Ball curve sðyÞ ¼ P 3 i¼0 P i S i ðyÞ have the following distinctive advantages over Bezier basis curve • The cubic Ball curve reduce to quadratic curve if P 1 = P 2 . • Operation of degree elevation and reduction for a polynomial curve with Ball basis functions are simpler and faster than the curve with Bezier basis functions [25].

GC 1 Rational Cubic Ball Curve
The cubic rational Ball curve s(θ) defined in [26,27] is where PðyÞ ¼ Að1 À yÞ 2 þ Bð1 À yÞ 2 y þ Cð1 À yÞy 2 þ Dy 2 ; QðyÞ ¼ ð1 À yÞ 2 þ að1 À yÞ 2 y þ bð1 À yÞy 2 þ y 2 : For GC 1 continuity consider the two consecutive curve segments where satisfying the following conditions using the above geometric continuity condition we can write where f i , f i + 1 , f i + 2 i = 1, . . ., n − 1 are the on curve control points and a i , b i and λ i are free parameters. For curve fitting, firstly, we will find off curve control points using least square method. Then, free parameters will be optimized with the help of genetic algorithm. In [28] the continuity between two curve segments is parametric continuity C 1 while in current work the continuity between two curve segments is geometric continuity GC 1 as the geometric continuity is more flexible than parametric. We have used C 1 rational Ball curves with tangents at end points. The tangent vectors work as intermediate control points, while in current work we are using GC 1 rational Ball curves and the intermediate control points are evaluated using least square method. In the current work the number of free parameters is increased to three and the proposed method is more flexible due to its geometric continuity. The proposed method works well for both small and large fractured parts. In current work we have extended the 2D CT scan data into 3D format taking equidistant z-component to construct the 3D craniofacial fractured part.

Effect of Control Point and Free Parameters
Ball curve can be changed or controlled by twofold. Firstly, by changing the off curve control points. For example in Fig 5, by changing the positron of P 1 control point the curve bend toward P 1 . Similarly, with P 2 control point keeping the free parameters unchanged. Secondly, by changing the free parameter a and keeping b constant. For example in Fig 6, the curve bends toward P 1 . Likewise, by changing b the curve will bend toward P 2 .
The given Dicom data are in 2D form. To construct the surface, patch one have to convert the data in 3D form. For this, we took the non decreasing z component of each 2D contours as  a height like z 0 < z 1 < . . . < z max . The i th contour is defined by the sequence of distinct data points which are counterclockwise ordered in the contour at the height z i .

Contours Blending
The fracture part curves for each 2D contour have been constructed using proposed interpolant, and is converted into 3D form taking equidistant z component [29]. Let C i − 1 (ψ), C i (ψ) and C i + 1 (ψ) are missing part curves at height z i − 1 , z i and z i + 1 respectively. The cubic Ball interpolant can be written as for i = 1, . . ., n − 1, 0 θ, ψ 1, where n represents the number of contours. S j , j = 0, . . ., 3. are Ball basis polynomial functions defined in Eq (1).
If λμ > 0, then Dividing denominator and numerator by ( This definition for v i (ψ) is proposed by [30] for a comonotone univariate interpolation scheme. Finally the surface between two consecutive contour z i and z i + 1 is given by, Now the surface between z 1 and z n is Sðc; yÞ ¼ fS i ðc; yÞ; i ¼ 1; :::; n À 1; 0⩽c; y⩽1g ð 10Þ

Boundary Extraction and Corner Detection
To construct the parietal bone fracture, we will initially determined the boundary of each CT scan image using mathematical morphology. The mathematical morphology is defined as β(A) = A − (AΘB). In this equation A represents the set of all black pixels, B represents 3 × 3 structured elements, β(A) is the boundary set of A, − and Θ represents the difference and erosion operator. To divide the boundary in smaller segments, we use the corner points. Sarfarz et al method [31] is employed to find the corner points.

Parameterization
This manuscript employs the chord length parametrization to find the values of θ i corresponding to D i , where D i represents the data points of each segments.
Normalized Mean Squares Error D i are the data points of each segments and θ is parameterized by chord length. We will find the free parameters a i , b i and λ i in proposed interpolant by minimizing the cost function. Genetic Algorithm (GA) defined in [28] is used to find the a i , b i and λ i respectively.

Graphical User Interface(GUI)
GUI consists of one or more windows having control, these are called components. GUI is utilized to perform the interactive tasks and display it graphically. GUI facilitates the users in completion of different tasks. Users are not required to learn the basic programming of each component. GUI consist of start and stop button, panel, scroll bar, push button and boxes etc.
GUI with its individual controls has a executable MATLAB code known as callbacks. The implementation of each callback is activated by a particular user action such as clicking a mouse button, pressing a screen button, typing a string or a numeric value, selecting a menu item, or passing the cursor over to a component. The GUI then reacts to these events. In this manuscript, GUI is used to construct and control the curves of the fractured segment of the parietal bone. The input for curve construction is missing part end points.

Proposed Algorithm
This section simplifies the algorithm for craniofacial fracture reconstruction.
Input: CT scan image in Dicom format. Output: Constructed 3D craniofacial fracture.  8. Conversion of 2D CT scan data into 3D form as in Fig 12   9. Construction of 3D craniofacial fracture using proposed method as in Figs 13 and 14.
In step 4 we reconstruct the boundary curve of a full skull using our proposed method. The genetic algorithm (GA) is used to optimize the free parameters. Normalized mean squares error is used as a cost function in GA and will be repeated until we get the best possible values of free parameters for boundary curve construction as shown in Fig 8(d). The process of optimization terminates when the error is less or equals to 10 −2 for the desired solution. In Fig 8  (d), reconstructed boundary curve is represented in red, while the original boundary curve is in black.

Case study: 3D Craniofacial Fracture Reconstruction
This section illustrates an example of application of the proposed algorithm. The given 2D CT scan Dicom data are in slices as in Fig 7. First, we constructed the boundary curves of complete skull using rational cubic Ball with GC 1 continuity at each knot as shown in Fig 8. Then, we have constructed the inner and outer curves of fractured part for all CT scan slices. Fig 9(a) is the original CT scan of slice 156. For the construction of fractured part, initially, we converted the original image to binary form as shown in Fig 9(b), then, using mathematical morphology, we extracted the boundary of skull as shown in Fig 9(c). Fig 9(d) represents the reconstructed inner and outer curves of fractured part using rational Ball curves. A least square method has been used to evaluate the intermediate control points of a rational Ball curve. The first and last control points of a rational Ball curve will be the starting and end point of the fractured part. The curves can be changed or alter using free parameters defined in Eq (6). Fig 10 shows the construction of fractured part curves of different slices. Fig 11 is the display of Graphical user interface (GUI) which helps in constructing the fractured curves. Next step is to convert the CT scan data and reconstructed fracture curves in 3D form taking equidistant z component as shown in Fig 12. Fig 12(a)

Conclusion
Contours blending function have been used to construct the 3D surface. A patient with parietal bone fracture have been discussed as a case study to show its applicability. Fracture part curves of each contour is constructed in 2D form using GC 1 rational Ball curves, since the given Dicom data are in 2D form. Then we convert each 2D contour into 3D form taking equidistant z component. The free parameters in proposed interpolant are optimized using Genetic Algorithm. The proposed method will provide the custom made implant for every individual patient. It is also time saving as the surgeon can perform the task of designing the implant geometry without the help of technicians and any other external help. It also reduces the incidence of infection. Proposed method is user friendly due to the presence of GUI.