Patient Specific Bone Modeling for Minimum Invasive Spine Surgery

This paper deals with patient specific bone modeling from non-invasive medical images viz. CT/MRI scans of scoliosis patients. The volumetric data i.e. voxel and triangular facet triangle based models are primarily used for bio-modeling and visualization, which requires huge memory space. On the other side, recent advances in CAD technology facilitate design, prototyping and manufacturing of any object having freeform surfaces. These CAD-based solid modeling is based on boundary representation (B-rep) techniques. Image Processing techniques are proposed to extract point cloud data from stalk of CT scan images. Estimated Point Cloud data is preprocessed and further used for water tight surface model. These patient specific bone models are more accurate, reliable and editable for implant design, trajectory path planning, etc. in robotic surgeries.


Introduction
Three-dimensional (3D) Bio-CAD model reconstruction from CT medical image has recently become the issue of much attention. This method is gaining popularity in Bio-medical engineering due to its capacity to create anatomical models with the help of technologies viz. Rapid Prototyping (RP) and Reverse Engineering (RE). These models have diagnostic, therapeutic and rehabilitatory medical applications. Bio-CAD is widely used in many applications such as Computer-Aided Surgery, Structural Modelling of Tissue, Design of Orthopaedic Device, Implants, Tissue Scaffolds and Freeform Fabrication or Bio-Manufacturing [1][2][3][4][5]. These models also have non-medical applications in the field of passenger safety design and crash analysis [3,[5][6][7][8]. The RE procedure consists of measure steps viz. Point Cloud Data capture, pre-processing of Point Cloud data (mainly for noise filtering), segmentation and CAD model creation [9,10]. Medical imaging modalities like Computer Tomography (CT), Magnetic Resonance Imaging (MRI) and Ultrasound are used for medical diagnostics and surgical planning. With the recent advances in computer graphics and scanning technologies, it is possible to generate 3D models for visualization or modelling [11,12]. The reconstructed 3D models can provide valuable medical information and powerful diagnostic tool for surgeons to understand the complex internal anatomy of the patient. These can also be used to fabricate prostheses and to perform various simulation and analytical tasks [2,9,11,12].
In 3D reconstruction of medical images, reconstruction methods are broadly classified into two major categories viz. the volume based approach and contour based approach. In majority of commercially available Bio-medical software, to interpret CT image data, volume pixels i.e. voxel representation is most commonly used. This volumetric model is represented by brick-like components, each representing a set height width, and depth. The voxel data from medical image data can be directly transferred to STL (stereo-lithography) format and printed on a rapid prototyping system. It is attractive in visualizing human anatomies since it can it generates high resolution surfaces and efficient in terms of computation and data structure. Voxel based representation suffers from problems like aliasing, loss of geometric information, huge memory requirements and inconvenient model editing. Also, triangular patched model from the marching cube is difficult for manufacturing and engineering applications without pre-processing [8,10].
These disadvantages of voxel models can be overcome by converting the CT images model into a CAD-based solid model. CAD-based solid modeling relies upon "boundary representation" techniques. In this method a solid model is defined by the boundaries that enclose it. These bounding surfaces are mathematically described by polynomial functions such as B-Spline curves, Non-Uniform Rational B-Spline (NURBS) functions, etc. It offers additional advantage of achieving Tangential (C 1 ) and Curvature (C 2 ) continuity over STL data, which offers only Positional (C 0 ) continuity. The boundary of an object is defined for every image in this approach either by tracing of the object or by edge detection by using image processing algorithm such as Canny, Previtt or Sobel Operator [9]. These CAD based methods offers the construction of the model by minimizing the size of the files and ensuring the closure of bounding surfaces. In addition to being a closer approximation to most 3D structures, the boundary represented CAD model is capable of being altered through Boolean operations andanalyzed by enabling CAD and analysis software packages.
In this paper methodology for development of patient specific implant design is discussed from non-invasive techniques.

Estimation of Point Cloud Data
Segmentation is the process in image analysis in which the object of interest is isolated from the background. The ultimate goal of segmentation is to identify the part of the data array that makes up an object in the real world. The role of segmentation discussed in this paper is to separate the bone of interest from its surroundings in a set of clinical CT images. The results can be used to identify regions containing the 3D surface of the bone for subsequent registration and image-based computer assisted orthopedic surgery [13].
Most of the current research work involves in collecting the feature points from images using spatial domain methods. Such information includes edges, interest or corner points, ridges etc. Segmentation is the process in image analysis in which the object of interest is isolated from the background. The ultimate goal of segmentation is to identify the part of the data array that makes up an object in the real world. Thresholding based techniques have been widely used in image segmentation [4,5,10,[13][14][15][16][17][18]. Image thresholding techniques can be classified into two categories viz. global thresholding and local thresholding. In global thresholding using manually or automatically determined threshold values, a point can be classified as object or background depending on its grey value. Since bone structures are of high intensity levels in CT images, they can usually be separated from soft tissue using thresholding based methods. Hangartner did a study to search an optimal threshold to assess bone parameters [10]. The partial volume to the partial volume effect, beam hardening due to polychromaticity of the X-ray beam, intensity in homogeneity of bone structures, and high grey level of surrounding pixels.
Local thresholding methods, on the other hand, calculate a different threshold value for each pixel based on local statistics. Some methods employ the mean plus standard deviation or the mean of the maximum and minimum values [4,5], median of pixel intensities [10,14], statistics based on local intensity gradient magnitude [10], but these methods may not perform well because of the partial volume effect and intensity in homogeneity. Li et al. presented a segmentation method using Bayesian relaxation labeling [3]. Wang et al. used the simplified Kang Method to validate bone segmentation with optional manual interaction to refine the resulting periosteal surface [19][20][21]. In these methods, the threshold values are chosen empirically. Maneka R. et al proposed use of Discrete Curvelet transformation and Features from Accelerated Segmentation Test (FAST) [4,5]. These approaches rely on initial thresholding to extract only the bony areas from each slice required for feature point extraction. Further Curvelet Transformation is applied to detect the edges from threshold slices. This paper proposes a fully automatic thresholding method for robust and accurate segmentation of the bone in CT images. First an initial segmentation is performed to divide the images into bone and non-bone classes. In the subsequent iterative process, by using image processing techniques like edge detection, boundaries of bones are identified. Figure 1 represents methodology adopted for generation of Point Cloud Data from CT scan images in DICOM 3.0 format by using automatic thresholding technique. A novel programming module is developed in MATLAB to get these point cloud data from input DICOM images.

Case study
The primary source of data for this work is CT images in DICOM 3.0 format as shown in Table 1. These grayscale images had been acquired from GE ProSpeed CT scanner in DICOM 3.0 format of 512 X 512 pixel size with data collection diameter 250 mm.
Considering upper left corner of image as origin, point cloud data is generated. Figure 2 shows stepwise results obtained from developed image processing module. From single CT scan slice, 2121 points are extracted. These points are stored as (x, y, z) coordinates in neutral file format 'ASCII' which can be further used in any CAD environment. This case study represents the first step of wider experimental protocol of conversion of CT images into 3D CAD model, which will provide higher order continuity (C 2 ) as well as provide additional advantages of CAD models like editing, etc. A new threshold based algorithm is designed for extraction of Point Cloud Data from CT scan. Since, bone is having heterogeneous structure, with variable density; estimation of threshold value is most important step in this algorithm. The proposed     use of these point cloud data is to construct B-Spline curves and surfaces which will enable us to achieve Positional (C 1 ) and Curvature (C 2 ) continuity. For rapid prototyping, this point cloud data can be transferred to *.stl file with only positional continuity (C 0 ). The major challenge in this work is huge number of points, to distinguish between inner and outer bone boundary profile points and detection of nearest points on one surface.

B-spline Curve Fitting
Various surface reconstruction algorithms like quadratic surface fitting, B-spline surface fitting, lofted surface fitting and sweep surface fitting methods have been proposed in the literature to fit 3D measured point cloud data into various surfaces [12,14,18]. To achieve Curvature Continuity (C 2 ), minimum three degree curve must be fitted through point cloud. Generally, for having three degree curve, Cubic splines, B-spline and NURBs can be used. Among them, B-spline surfaces, especially Non-Uniform Rational B-Spline (NURBs) surfaces, are popular due to their ability to accurately approximate most types of surface entities encountered in design and manufacturing application [14]. Dong-Jin Yoo proposed a method to reconstruct a B-spline surface from a point cloud data or a sequence of CT image data is proposed by using an efficient implicit surface interpolation scheme [14]. Closed Non-Uniform Rational B-Spline (NURBS) curves can be used to represent the contours of the layers to maintain the surface accuracy of the CAD model [3]. Olya Grove et al proposed use of medical image data is to reconstruct a surface from serial parallel contours extracted from images [18]. Data points representing the contours of the region interest are extracted and fitted with B-Spline curves, and the sequence of curves is lofted to generate the 3D surface model. These control point based of methods viz. B-splines and NURBs provides excellent coverage but high degree of freedom due to large number of control points. These curves provide higher order continuities as compared to "stl" and voxel based methods.

Preprocessing of point cloud data
CT scan is a stalk of equispaced images, the generated points are in equispaced layers of 5 mm. Figure 3 shows a dense point cloud data generated from DICOM image. This data contains points of different edges viz. external boundaries, internal boundaries, intra bone cavities, etc. Image processing algorithm uses line scan approach for storage of point cloud data.
For automatic curve generation, following issues must be tackled for fitting curves through generated Point Cloud Data -For single loop passing through points, all points are needed to be arranged in a sequence of nearest point. Figure 4 shows implementation of Euclidian nearest point algorithm developed for detection of nearest point in a loop and arranging them into a proper sequence.
For two loops i.e. internal and external boundaries, if curve fitting implemented, it results into zig-zag pattern. A Min-max algorithm has been developed to separate inner and outer points as shown in Figure 5. To overcome this issue, a tracking algorithm has been developed as shown in Figure 7 for sorting points for each loop separately.

Tracking Algorithm
Convert CT images into binary images.     Then Scan for value equal to threshold value in its 8 neighbours. There will be 2 positions where its value is 1. So, select any one in order to have curve clockwise or anti-clockwise.
Change the value of S (i, j) = 0. So, that point is not repeated or reiterated.
Again iterate with new (i, j) value & store its position in matrix T.
Once loop is over, store value in new matrix with new curve's pixel position.
Thus, organized points are achieved & pass spline curve through each matrix points.

Fitting B-Spline Curve
Conventionally, interpolation and approximation [22] are the two widely used techniques used to fit a curve to a sizable data set. In interpolation, i.e. the curve passes through all the data points, and then removes as many knots (control points) as allowed by the tolerance. This method did not work because the interpolating curve had too many wiggles that were not removable even with a high tolerance. In second approach, an approximation with some number of control points and increase the degree of freedom (knots and control points) until the tolerance requirement is met. It fails due to failure of estimating degree of freedom, also it loses intricate shapes. To improve the modelling process, Point Cloud Data smoothening is required before processing into CAD environment. Figure 8 shows a sample Point Cloud Data and closed B-spline curve fitting through sample point cloud data by using Imageware software. Figure 9 shows process of conversion of CT scan of 5 years female suffering from diastematomyelia resulting into splitting spinal cord in two parts in spinal column to rapid prototyped model. Figure 10 shows RP models developed from CT scan of various scoliosis patients.

Conclusion
This paper discussed initial steps of wider experimental protocol of conversion of CT images into 3D CAD model, which will provide higher order continuity (C 2 ) as well as provide additional advantages of CAD models like editing, etc. A new threshold based algorithm is designed for extraction of Point Cloud Data from CT scan. By using Image Processing techniques for noise reduction and edge detection, edges were identified, which is further converted into (x, y, z) coordinates and data stored in ASCII file, which can be further used in any CAD environment. The major challenge in this work is huge number of points, to distinguish between inner and outer bone boundary profile points, detection of nearest points on one surface. NURB's surface with Curvature (C 2 ) continuity is fitted through these point cloud data.