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].

Fig 1

Craniofacial bones.

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-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 (GC1) 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.

Cubic Ball Basis Functions and Curve

The cubic Ball polynomial basis was first proposed by Ball [24] for CAD systems application. Fig 2 illustrates these functions against its parameter θ. The Ball basis functions can be written as

Fig 2

Ball basis functions.

Theorem 1: The Ball basis functions defined in Eq (1) have the following properties

Linearly Independent: The Ball basis functions are linearly independent. There do not exist set of non zero constants a0, …, a3, for which.

Non Negative: Ball basis functions are non negative when θ ∈ [0, 1].

Symmetric: The Ball basis functions are symmetric as

Monotonicity: S0(θ) is monotonically decreasing and S3(θ) is monotonically increasing when θ ∈ [0, 1].

Partition of Unity: Ball basis functions form a partition of unity, it means that sum of Ball basis functions will be 1.

Theorem 2: Let P be the set of control points and S(θ), i = 0, …, 3 are Ball basis functions defined in Eq (1). The Ball curve 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.

Variation Diminishing Property (VDP): Variation Dimension Property is obeyed by Ball curves. VDP is stated as if a curve is intercepted by straight line in b number of points and control polygon by p number of points, then it will always hold that as shown in Fig 3. Where k is 0 or a positive integer.

Fig 3

Variation diminishing property with b = 1 p = 3 k = 1.

Endpoint Interpolation: The Ball curves always passes through the first and last control points.

Theorem 3: Let P be the set of control points and S(θ), i = 0, …, 3 are Ball basis functions defined in Eq (1). The Ball curve have the following distinctive advantages over Bezier basis curve

The cubic Ball curve reduce to quadratic curve if P1 = P2.

The Ball curve passed away from the convex hull intermediate control points as compared to the Bezier curve as shown in Fig 4.

Fig 4

Ball verses Bezier curve.

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].

GC1 Rational Cubic Ball Curve

The cubic rational Ball curve s(θ) defined in [26, 27] is where For GC1 continuity consider the two consecutive curve segments where satisfying the following conditions using the above geometric continuity condition we can write where f, f, f i = 1, …, n − 1 are the on curve control points and a, b and λ 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 C1 while in current work the continuity between two curve segments is geometric continuity GC1 as the geometric continuity is more flexible than parametric. We have used C1 rational Ball curves with tangents at end points. The tangent vectors work as intermediate control points, while in current work we are using GC1 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 P1 control point the curve bend toward P1. Similarly, with P2 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 P1. Likewise, by changing b the curve will bend toward P2.

Fig 5

Effect of control points on curve.

Fig 6

Effect of free parameters black curve when a = 1, b = 1 blue a = 1.4, b = 1 black a = 1.6, b = 1 cyan a = 1.8, b = 1 red a = 2, b = 1 black.

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 z0 < z1 < … < z. The i contour is defined by the sequence of distinct data points which are counterclockwise ordered in the contour at the height z.

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(ψ), C(ψ) and C(ψ) are missing part curves at height z, z and z 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 = 0, …, 3. are Ball basis polynomial functions defined in Eq (1).

v(ψ) is a tangent to surface at (C(ψ), z), to define the tangent vector let

P(ψ) = (C(ψ), z) for i = 1, …, n. Where P(ψ) is in 3D form because C(ψ) is in 2D form. Let Thus is a complex combination of the vector (P(ψ) − P(ψ)) and (P(ψ) − P(ψ)). The z − component of is

|C(ψ) − C(ψ)|(z − z) + |C(ψ) − C(ψ)|(z − z). Let V(ψ) be the vector product when is divided by its z − component.

Then V(ψ) = (v(ψ), 1) where where D1 = C(ψ) − C(ψ) and D2 = C(ψ) − C(ψ).

The vector V(ψ) will be chosen to be a tangent to the surface at P(ψ). If C(ψ) = C(ψ) = C(ψ), then D1 = D2 = 0 and v(ψ) is set to be (0, 0).

As v(ψ) lies in a plane containing P(ψ), P(ψ) and P(ψ). If these points lie in a vertical plane, then for some unit vector v*(ψ) and some scalers λ and μ. Then If λμ < 0, then |μ|λ = −|λ|μ and v(ψ) = 0.

If λμ > 0, then Dividing denominator and numerator by (z − z)(z − z), This definition for v(ψ) is proposed by [30] for a comonotone univariate interpolation scheme.

Finally the surface between two consecutive contour z and z is given by, Now the surface between z1 and z is

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 θ corresponding to D, where D represents the data points of each segments.

Normalized Mean Squares Error

D are the data points of each segments and θ is parameterized by chord length. We will find the free parameters a, b and λ in proposed interpolant by minimizing the cost function. Genetic Algorithm (GA) defined in [28] is used to find the a, b and λ 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.

Read image as in Fig 7.

Fig 7

Given CT scanned images of patient with parietal bone fracture.

Boundary extraction Fig 8(c)

Fig 8

Reconstructed image of CT scan of Dicom slice 178.

Corner points detection to make segments of curve Fig 8(d)

GC1 rational Ball curve is used to fit each segment. The unknown parameters a, b and λ in Eq (2) are optimized using Genetic Algorithm Fig 8(d)

Step 4 is repeated until a desired solution is obtained

Fracture part curves reconstruction for each CT scan slice as in Figs 9 and 10

Fig 9

CT scanned data of slice 156.

Fig 10

Reconstructed missing part curves for different CT scanned slices.

Using GUI for the construction of curves as in Fig 11

Fig 11

GUI display of craniofacial reconstruction.

Conversion of 2D CT scan data into 3D form as in Fig 12

Fig 12

CT scan data in 3D format.

Construction of 3D craniofacial fracture using proposed method as in Figs 13 and 14.

Fig 13

The front view of reconstructed 3D craniofacial fracture.

Fig 14

The side view of reconstructed 3D craniofacial fracture.

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 GC1 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) and 12(b) are front and side view of CT scan data in 3D form respectively, similarly Fig 12(c) and 12(d) represent the front and side view of 3D CT scan data with fractured part curves respectively. Figs 13 and 14 represent the front and side view of constructed 3D parietal bone fracture using proposed contour blending function.

The Human Research Ethics Committee, Universiti Sains Malaysia, approved the study. Data were collected from archived images of patients from the Picture Archiving and Communication System (PACS) server at the Radiology Department, Hospital Universiti Sains Malaysia. The patients had undergone decompressive craniectomy surgery between May 2009 and December 2014 at the Hospital Universiti Sains Malaysia, Kubang Kerian, Kelantan Malaysia. The medical images were anonymized and the researchers had no access to patients’ information.

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 GC1 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.