Abstract
The Sensimmer platform represents our ongoing research on simultaneous haptics and graphics rendering of 3D models. For simulation of medical and surgical procedures using Sensimmer, 3D models must be obtained from medical imaging data, such as magnetic resonance imaging (MRI) or computed tomography (CT). Image segmentation techniques are used to determine the anatomies of interest from the images. 3D models are obtained from segmentation and their triangle reduction is required for graphics and haptics rendering. This paper focuses on creating 3D models by automating the segmentation of CT images based on the pixel contrast for integrating the interface between Sensimmer and medical imaging devices, using the volumetric approach, Hough transform method, and manual centering method. Hence, automating the process has reduced the segmentation time by 56.35% while maintaining the same accuracy of the output at ±2 voxels.
Keywords: Image segmentation, Haptic rendering, Surgical simulation, Computed tomography (CT), DICOM
Introduction
Model-based segmentation in medical imagery has undergone rapid advancements over the past decade. Computer tomography (CT) is the most common data type used in medical imagery. The image segmentation helps in generating labels for each set of pixels in the image; hence, by using a stack of CT images, a 3D model can be reconstructed. There is a complexity in handling noisy datasets to differentiate the region of interest from the noise in the image.
Sensimmer represents the ongoing evolution of ImmersiveTouch [1, 2], the latest generation of augmented virtual reality (VR) technology, which integrates a haptic device with a head and hand tracking system and a high-resolution, high-pixel-density stereoscopic display. A haptic device collocated with 3D graphics is the key to deliver extremely realistic simulations. Previously, ImmersiveTouch has been successfully applied to the simulation of neurosurgical procedures and training of neurosurgery residents [3]. It implements graphics and haptics rendering in a multithreaded environment.
To satisfy the minimum required graphic and haptic frame rates, it is essential to use efficient 3D models of the anatomical region to be simulated. This efficiency involves a trade-off. It is necessary to maintain an adequate polygon count in order not to lose any significant detail, but polygon reduction with re-topology must also be performed to ensure the desired haptic frame rate. In the context of on-demand high-fidelity simulations [4], automatic or semi-automatic techniques to generate 3D models from medical images are highly desirable. Segmentation techniques are applied to magnetic resonance imaging (MRI) or computed tomography (CT) images to obtain 3D models from medical data. These models must be reduced and converted to polygonal surfaces for simultaneous graphics and haptics rendering. Depending on the application, further processing may be needed (e.g., drilling of burr holes for simulation of a neurosurgery procedure known as ventriculostomy).
This paper describes ongoing research on the process of generating the 3D model to be used in Sensimmer. For this purpose, many freely available open-source software tools are evaluated, namely, the Insight Segmentation and Registration Toolkit (ITK) [5], ITK-SNAP [6], MATLAB [7], and Seg3D [8]. Using these tools, the segmentation of various anatomical features such as the vertebrae, discs, spinal cord, and the aorta can be done. Since manual segmentation is time consuming, a new semi-automated segmentation method has been developed. The results are evaluated with the goal of modeling 3D anatomy to within ±2 voxels of the original image. Moreover, the choice of method depends on the grayscale values and visibility of the anatomies of the CT dataset.
This paper is organized as follows. “Introduction” provides a brief review of a few papers selected for their use of particular methods of interest. “Background” describes the detection of the boundary and shape for the desired anatomy with the generalized Hough transform (GHT). “Methods” describes the methodology of segmentation based on identifying the different contrast values that correspond to the various anatomical structures regardless of the source of the image. “Process Evaluation” shows the comparison of the result based on the time considering the same accuracy. The conclusion is presented in “Conclusion.”
Background
A detailed literature survey was performed to locate various works on medical image segmentation. However, since our demonstration concerns automatic image segmentation of anatomical features that resemble an ellipse or circle, such as the vertebrae, intervertebral discs, spinal cord, and the aorta, we have limited our report here to papers that deal with segmentation of the discs and the aorta.
Xie et al. [9] explored the concept of detection of elliptical shapes. The method of detecting a circular shape was considered as the starting point, and the capability of elliptical shape detection was developed using a modified Hough transform [10]. Primarily three parameters were considered for detecting a circular shaped object, i.e., radius, circumference, and the center point. A similar method was adopted for an ellipse but considers five parameters, i.e., the circumference, major and minor axis, and two points which are initially considered on the ellipse, to obtain the loop around which the ellipse is formed. Adapting this paper for the present study, circles were used to detect the spinal cord and aorta, whereas an ellipse was used to detect the intervertebral discs.
Neubert et al. [11] proposed an automated approach to extract the 3D model of lumbar and thoracic intravertebral discs (IVDs) and vertebral bodies (VBs) from MR images using statistical shape analysis and registration of gray level intensity profiles [12]. The algorithm was divided into two main parts: (1) spine localization and (2) shape-based segmentation of IVDs and VBs. During the spine localization, the 3D spine curve was extracted and approximate positions of VBs were determined using active rectangles. The shapes were subsequently deformed using the active shape modeling (ASM) strategy [13].
Behrens et al. [14] illustrated the use of randomized Hough transform (RHT). They presented a new approach for the coarse segmentation of tubular structures in 3D image data. The algorithm requires only a few initial values and minimal user intervention and can be used to initialize complex deformable models. It is based on an extension of the RHT. This paper helped us to obtain an insight on how to extrude the elliptical 2D image into a 3D cylindrical structure. We adopted the concept of detecting the elliptical structure from a 5-point coplanar method and used a vector equation of the ellipse to obtain the cylinder axis.
Wang et al. [15] present an automatic method to segment spinal canals in low-resolution, low-contrast CT images. The scans were collected from eight different sites and varied significantly in field-of-view (FOV), resolution, pathology, etc. This paper adapted the idea of interactive segmentation to form a fully automatic approach that segments spinal canals from CT images. To start the automatic pipeline, the authors identified voxels that are inside the spinal canal according to their appearance features. Then the detected seeds were input to random walks (RW) to produce the segmentation of foreground/background. Seeds were adjusted accordingly and fed back to RW for better segmentation. RW asks users to specify seeding voxels of different labels, and then assigns labels to non-seeding voxels by embedding the image into a graph and utilizing intensity similarity between voxels. Users can edit the placement of seeds to acquire more satisfactory results. Based on analysis, iterative optimization successfully enhanced the capability of RW in dealing with tubular spinal canals and showed that the boundary conditions can be improved to guarantee better segmentation results.
Methods
The insight segmentation and registration toolkit (ITK) [6, 16] is an application framework initially developed to support the U.S. National Library of Medicine’s Visible Human Project [17]. ITK-SNAP [6] is an open-source software package, built on top of ITK, oriented to the segmentation of 3D anatomical structures from medical images. Using ITK-SNAP, it is possible to perform segmentation as a semi-automated procedure. Though referred to as snakes within the software, ITK-SNAP uses two 3D active contour segmentation methods [18], namely, geodesic active contours [19], driven by intensity edges, and region competition [20], driven by intensity regions. Similar tools including MATLAB and Seg3D were later identified and proved to reduce the amount of time in segmenting the required anatomy compared to ITK-SNAP.
MATLAB is a programming platform with a proprietary programming language developed by MathWorks. The image processing toolbox in MATLAB allows image processing and analysis by using the available algorithms and functions and modifying them for specific cases. Seg3D is a volume segmentation and processing tool created by the NIH NCRR Center for Integrative Biomedical Computing (CIBC) located at the Scientific Computing and Imaging Institute (SCI) at the University of Utah, and was developed in collaboration with Numira Biosciences. It combines a flexible manual segmentation interface with powerful higher dimensional image processing and segmentation algorithms from the Insight Toolkit. Seg3D data is handled mainly in the form of projects which are similar to those used in other applications. A project consists of a group of files, and with this setup, Seg3D can track and save data that we are working on, the settings that we are using, and the tools and filters.
Below, we present the methods used to extract different models from four different datasets. The segmentation of the anatomy is categorized into two types, the volume model and the surface model. The vertebrae and disc fall under the volume models and the spine and aorta are surface models, based on the anatomy of the organ. Each dataset is unique in terms of the number of CT images, the anatomy, and its visibility.
To evaluate the segmentation accuracy based on the anatomy, object recognition is done to identify the shape and detect the boundary [21]. Hence, we use the shape detection methods to understand the accuracy of the segmentation of the CT images.
Object Detection
The GHT [22, 23] can be used to detect the boundaries of the desired anatomy in an image. To achieve our results, we used the Hough transform feature within MATLAB; but in this section, we introduce the GHT as such. The GHT is an extension of the Hough transform (HT) [23, 24], a standard template-matching algorithm. In general, an anatomical area of interest is represented as an array of landmark points within the area of interest, H(x i,y i), of the same size as the input image. The array consists of a 2D accumulator to represent the CT slice. Each pixel value of H(x i,y i) determines whether a given pixel lies within the area of interest according to predetermined grayscale values that define the threshold. The anatomical reconstruction is initiated by an operator who chooses an arbitrary seeding point within the anatomical feature (e.g., disc, aorta), which allows the full shape to be grown by iteration: each pixel edge in the image “votes” for the next pixel to be either included in or excluded from the region of interest. (It should be noted that the aorta is grown as a solid volume in our work.) To perform the GHT, the template shapes are built in advance, requiring the user to know exactly what shapes will be encountered. Further, when the scale and orientation of an input shape are variant and unknown in advance, brute force is usually employed to enumerate all possible scales (S i) and orientations (Θ i) of the input shape in the GHT process. This adds two dimensions to the parameter space, thus requiring a 4D accumulator, H(x i,y i,S i,Θ i). This dramatically increases the execution time and leads to sparsity in the accumulator, making the selection of strong matches more difficult. Iterative approaches to the GHT, as proposed in [25], eliminate the extra dimensions by using a local accumulator to find each template in turn and a global accumulator to collect the best score from the local accumulator.
Shape Detection
In order to model a scale- and orientation-invariant shape, we represent it by a set of points (r i,θ i) in polar coordinates, as in Ref. [25]:
Assume that OA is an arbitrary radius of the template. In our algorithm, starting from OA and moving clockwise, we divide the circle into n equal arcs to place points around the boundary, with each arc being (360/n) degrees. So, the above set of points can be represented as
and the corresponding vector p in x-y coordinates consists of the points
Segmentation Strategy: MATLAB and Seg3D
A key finding of this research is that the methods we have devised in MATLAB produce better segmentation outcomes than those from Seg3D in those cases where the datasets include a representation of a given anatomical feature that is without noise, e.g., a complete cross section of a disc or the aorta. Here, MATLAB has the advantage because it can complete the segmentation without manual intervention. By contrast, Seg3D is a semi-automatic approach. The user visually inspects the image and then places “seeds” on the image to serve as starting points from which the algorithm can find the pixels with comparable contrast values to generate the entire anatomical feature, e.g., the entire aorta or one or several intervertebral discs.
Visual inspection by a trained operator determines whether the visibility of a given anatomical feature is “good” and can be reconstructed by MATLAB, or “average” or “bad” and should be processed using Seg3D. Accordingly, in Table 1 below, datasets 1–4 for the vertebrae, 3 for discs, and 1–2 for the aorta were processed by MATLAB, and the rest by Seg3D. Table 2 shows the contrast level for the best possible visibility of each respective anatomical feature. Every image is in DICOM format and contains intensity levels of 512 by 512 pixels. Intensity levels are given as 16-bit signed integers.
Table 1.
Dataset # | Vertebrae | Discs | Spinal cord | Aorta |
---|---|---|---|---|
Dataset 1 | Good | Average | Average | Good |
Dataset 2 | Good | Average | Bad | Good |
Dataset 3 | Good | Good | Bad | Bad |
Dataset 4 | Average | Bad | Bad | Bad |
Table 2.
Dataset # | Vertebrae | Disc | Spinal cord | Aorta |
---|---|---|---|---|
Dataset 1 | 1214 ± 166 | 1050 ± 100 | 1200 ± 100 | 1139 ± 20 |
Dataset 2 | 1146 ± 200 | 1050 ± 50 | 1250 ± 50 | 1206 ± 13 |
Dataset 3 | 1142 ± 4 | 1075 ± 100 | 1405 ± 300 | 1013 ± 10 |
Dataset 4 | 1104 ± 49 | 1100 ± 25 | 1249 ± 100 | 976 ± 15 |
The contrast ranges of the different anatomical features depend on the CT machine used and the method followed during the scanning process. The contrast values tabulated above show the best possible visibility of the respective anatomy in the test datasets used and were obtained using MATLAB. The software has a function in its image processing toolbox which enables the user to adjust the histogram and vary the contrast of the image. This proved to be very useful in identifying the suitable contrast range for the best visibility of the anatomy.
The following figure shows the process flow for the datasets with good, average, and bad contrasts using MATLAB and Seg3D (Figs. 1 and 2).
Segmentation Using MATLAB
Segmentation using MATLAB starts with importing the DICOM files into MATLAB, organizing the files in sequence, and normalizing the images to the contrast range that provides the best possible visibility of the organ to be segmented. The approach varies based on these three different methods, namely, the volumetric approach, manual centering method, and HT method.
In the volumetric approach, the algorithm detects volumes and then segments them based on set contrast differences. Since the most common form of CT images has good contrast, the vertebrae appear the brightest. We apply the volumetric method only to the vertebrae, not to soft tissues. This method detects vertebrae along with other bones in the images. The other bones are masked to leave them black and thus exclude them from the region of interest containing the vertebrae. This approach gives an output (a volume) that is not governed by any functions, does not have surface points, and cannot be meshed.
A variation of the volumetric approach is the watershed algorithm. The watershed algorithm transforms the gradient of a gray level image into a topographic surface. The algorithm punctures holes at the local minimum of the intensity and fills the region with water. There are three stages in our application of the watershed algorithm. First, a directed graph searching algorithm is applied to find the longest path from source to sink, which is the spinal canal in our case. The vertebral column segmentation is implemented using geometric and statistical models owing to its articulated and complex structure. The last stage is to remove the spinal cord from within the spinal canal. We developed a partition approach based on the intensity profile of the reformatted images. The centerline of the spinal canal is first projected onto the reformatted images. Then the normal is computed at every point on the centerline. The intensity along the normal direction is aggregated and recorded. Therefore, the watershed algorithm can be used only for vertebrae and discs, but not for the spinal cord and the aorta, which must be surface models (Fig. 3).
The manual centering method is a semi-automatic method where the user needs to click on the center of the organ and a selected area around it is masked. It is applied here to the segmentation of the spinal cord and aorta. The image processing takes place within the masked region based on contrast differences. The organ is identified as its edges are detected, its interior points are obtained, and finally a mesh is created after the algorithm has been run on a loop, for every few slices, for the entire dataset. This method is highly suitable for spinal cord and aorta, since these two organs are almost circular and have a regular shape throughout the dataset. Difficulties identified in this method included the leakage in the mesh, which causes the final output to mutate outside the organ’s actual position, and the fact that the user needs to manually enter the center points of the organ once every few slices, which takes a long time for larger datasets.
Figs. 4 and 5 show the point selected and the subsequent result obtained from the manual centering method by detecting the circle around the marked center. The coordinates of centers and radii of all slices are obtained by interpolating the available coordinates, and then the 3D surface mesh is created.
Another approach available in MATLAB is the HT method, which is an algorithm that is run to detect certain shapes based on the parameters given by the user. This algorithm runs the detection code in a loop for the entire dataset. The method was used to detect circles for segmenting the spinal cord and/or aorta and ellipses for segmenting the discs. A variety of functions were used to tailor the general code for the specific anatomy of a given dataset. The result of this method is also an array containing the coordinates of the centers and the radii, which was plotted and used to create a surface mesh. Disadvantages of this method were that the algorithm failed to detect the shape in some cases, and the code had to be modified for each anatomy of each dataset, which is a time consuming activity. Also, the algorithm was not able to detect the shape in the branching region of the aorta. Fig. 6 shows the result obtained when segmenting a non-branching aorta of good contrast using the HT.
Segmentation Using Seg3D
The manual segmentation technique is used for identifying the anatomy visually and segmenting the area of interest. Here, the segmentation is done using Seg3D. The segmentation in Seg3D starts with importing the series of images into Seg3D as data volume. The images appear as a single data layer over which subsequent operations can be performed. Seg3D has two types of layers that can be operated by different tools and filters: (i) data layer and (ii) mask layer. Data layer is the image data that has been imported, or the result of Gaussian, mean, or median filters. The mask layer displays the segmented regions of interest as a colored opaque layer. Otsu threshold filter was used to display the image intensity histogram to determine adequate threshold levels and segment the image based on the determined threshold levels. It was necessary to use the threshold tool on the dataset to segment the desired anatomy. So, the data noise in the image was suppressed using a Gaussian blur filter. But the output from this thresholding filter was not completely accurate and required a lot of manual cleanup. For doing the manual cleanup, the polyline and paintbrush tools were used. The polyline tool enabled filling/removal of segmented areas at a faster rate compared to the paint brush tool, which was used only for making precision cleanup. The availability of logical operators in Seg3D proved to be highly useful in operations requiring a combination of two or more layers. Apart from this thresholding, the other method that was effective in segmenting the required parts was the confidence connected filter. In this method, the criterion is based on simple statistics of the current region. First, the algorithm computes the mean and standard deviation of intensity values for all the pixels currently included in the region. A user-provided factor is used to multiply the standard deviation and define a range around the mean. Neighbor pixels whose intensity values fall inside the range are accepted and included in the region. When no more neighbor pixels are found that satisfy the criterion, the algorithm is considered to have finished its first iteration. At that point, the mean and standard deviation of the intensity levels are recomputed using all the pixels currently included in the region. This mean and standard deviation define a new intensity range that is used to visit current region neighbors and evaluate whether their intensity falls inside the range. This iterative process is repeated until no more pixels are added or the maximum number of iterations is reached. Finally, dilate and erode filters were applied to the output of the previous step and this filter was very useful in smoothing a mask layer, as performing a dilation and erosion on a mask volume will fill in some of the details on the surface.
Integration of All Models into a Single File
The integration of all models into a single file was done using MATLAB. The outputs were obtained either as a pair of .mhd and .raw files or as DICOM files. These outputs were opened and combined into a single file using MATLAB. The integrated output is used for additional processing before being used in the simulator. The following figure shows the integrated output viewed in ITK-SNAP.
Process Evaluation
To evaluate the effectiveness of different process flows, the total time required for segmentation was considered with an accuracy of ±2 voxels. The resulting comparison of all dataset is shown in Figs. 7 and 8. Here, it can be noted that for each dataset, the time taken for manual and semi-automatic approaches differ; this is because the contrast intensity of the dataset for each organ differs as shown in Tables 1 and 2.
It is to be noted that although the semi-automatic methods we developed have many steps compared to manual segmentation, the overall time taken is less for good datasets. The following Gantt charts show the time taken for segmenting different anatomies in different methods. The chart (Fig. 9 and Fig. 10) shows that the manual segmentation pipeline for dataset 4 takes 2433 min in total, while the proposed approach we developed takes only 1062 min to finish segmenting the anatomies with the accuracy considered the same as the manual segmentation method.
Conclusion
Extraction of accurate and efficient 3D models from medical images is a fundamental part of our VR and haptics simulator. Due to the large variability in different patients, this modeling component is difficult to automate completely. However, a reasonable degree of automation has been achieved using MATLAB and Seg3D. This paper describes the segmentation of four datasets by semi-automatic approaches like the volumetric approach, manual centering method, and HT method, all of which proved their suitability to solve the specific problem of time consumption with the same degree of accuracy as manual segmentation. Nevertheless, the search for a fully automated segmentation process, for example with the statistical shape model, remains ongoing.
Acknowledgement
Dr. Daniel P. Bailey, manager of research planning, University of Illinois at Chicago College of Engineering, contributed substantially to the writing and copy-editing of this article.
References
- 1.Luciano C, et al: Design of the ImmersiveTouch™: a high-performance haptic augmented virtual reality system. 11th International Conference on Human-Computer Interaction, Las Vegas, NV, 2005
- 2.Banerjee P, et al: Compact haptic and augmented virtual reality system. U.S. Patent No. 7,812,815. 12 Oct. 2010
- 3.Luciano C, et al. Second generation haptic ventriculostomy simulator using the ImmersiveTouch™ system. Studies in health technology and informatics. 2005;119:343. [PubMed] [Google Scholar]
- 4.Banerjee P, Charbel F: On-demand high fidelity neurosurgical procedure simulator prototype at University of Illinois using virtual reality and haptics. Accreditation Council for Graduate Medical Education (ACGME) Bulletin 20–21, 2006
- 5.The visualization toolkit. Available: http://www.kitware.com
- 6.The insight segmentation and registration toolkit. Available: http://www.itk.org/
- 7.MATLAB support. Available: http://www.mathworks.com/help/matlab/
- 8.Seg3D. Available: http://www.sci.utah.edu/cibc-software/seg3d.html
- 9.Xie Y, Ji Q: A new efficient ellipse detection method. Pattern Recognition, 2002. Proceedings. 16th International Conference on. Vol. 2. IEEE, 2002
- 10.Singh S, Datar A. EDGE detection techniques using Hough transform. International Journal of Emerging Technology and Advanced Engineering. 2013;3.6:333–337. [Google Scholar]
- 11.Neubert A, et al. Automated detection, 3D segmentation and analysis of high resolution spine MR images using statistical shape models. Physics in medicine and biology. 2012;57.24:8357. doi: 10.1088/0031-9155/57/24/8357. [DOI] [PubMed] [Google Scholar]
- 12.Davies R, et al. 3D statistical shape models using direct optimisation of description length. Computer vision—ECCV. 2002;2002:1–17. [Google Scholar]
- 13.Heimann T, Meinzer H-P. Statistical shape models for 3D medical image segmentation: a review. Medical image analysis. 2009;13.4:543–563. doi: 10.1016/j.media.2009.05.004. [DOI] [PubMed] [Google Scholar]
- 14.Behrens T, Rohr K, and Siegfried Stiehl H: Using an extended Hough transform combined with a Kalman filter to segment tubular structures in 3D medical images. VMV 2001
- 15.Wang Q et al: Automatic and reliable segmentation of spinal canals in low-resolution, low-contrast CT images. Computational methods and clinical applications for spine imaging. Springer International Publishing 15–24, 2014
- 16.ITK-SNAP. Available: http://www.itksnap.org/
- 17.The visible human project. Available: http://www.nlm.nih.gov/research/visible/visible_human.html
- 18.Yushkevich P, et al: Geodesic snakes for user-guided segmentation of 3-D anatomical objects: significantly improved efficiency and reliability. Neuroimage 2005 [DOI] [PubMed]
- 19.Caselles V, Kimmel R, Sapiro G. Geodesic active contours. Int J Comput Vis. 1997;22.1:61–79. doi: 10.1023/A:1007979827043. [DOI] [Google Scholar]
- 20.Rizzi SH, Banerjee PP, and Luciano CJ: Automating the extraction of 3D models from medical images for virtual reality and haptic simulations. Automation Science and Engineering, 2007. CASE 2007. IEEE International Conference on. IEEE, 2007
- 21.Lee K-M, Street WN: Learning shapes for automatic image segmentation. Proc. INFORMS-KORMS Conference. 2000
- 22.Ballard DH. Generalizing the Hough transform to detect arbitrary shapes. Pattern recognition. 1981;13.2:111–122. doi: 10.1016/0031-3203(81)90009-1. [DOI] [Google Scholar]
- 23.Chung C-H, Cheng S-C, Chang C-C. Adaptive image segmentation for region-based object retrieval using generalized Hough transform. Pattern recognition. 2010;43.10:3219–3232. doi: 10.1016/j.patcog.2010.04.022. [DOI] [Google Scholar]
- 24.Hough Paul VC: Method and means for recognizing complex patterns. U.S. Patent No. 3,069,654. 1962
- 25.Lee K-M, Street WN, Lee K: A fast and robust approach for automated segmentation of breast cancer nuclei. In Proceedings of the IASTED International Conference on Computer Graphics and Imaging. 1999