A multi-scale method for automatically extracting the dominant features of cervical vertebrae in CT images

Localization of the dominant points of cervical spines in medical images is important for improving the medical automation in clinical head and neck applications. In order to automatically identify the dominant points of cervical vertebrae in neck CT images with precision, we propose a method based on multi-scale contour analysis to analyzing the deformable shape of spines. To extract the spine contour, we introduce a method to automatically generate the initial contour of the spine shape, and the distance field for level set active contour iterations can also be deduced. In the shape analysis stage, we at first coarsely segment the extracted contour with zero-crossing points of the curvature based on the analysis with curvature scale space, and the spine shape is modeled with the analysis of curvature scale space. Then, each segmented curve is analyzed geometrically based on the turning angle property at different scales, and the local extreme points are extracted and verified as the dominant feature points. The vertices of the shape contour are approximately derived with the analysis at coarse scale, and then adjusted precisely at fine scale. Consequently, the results of experiment show that we approach a success rate of 93.4% and accuracy of 0.37mm by comparing with the manual results.


I. INTRODUCTION
Anatomical landmarks and dominant points of cervical vertebrae are of considerable importance for many applications on orthopedics, neurology, and radiation therapy planning.In many researches about computer-aided techniques, the geometric characteristics of anatomical features were mentioned to be utilized on the applications like image-based surgical guidance and operation planning [1] .However, because of the anatomical variation between patients and the complexity of medical images, automatic analysis and information collection in computerized tomography (CT) images are still challenging tasks.
In [1] , Lee et al. proposed a method to automatically locate the lumbar spine pedicles in CT images by referencing the canal boundaries for pedicle screw.But for cervical vertebrae operation planning, more landmarks are required [7] .Dominant feature points of cervical vertebrae include transverse foramens, spinous processes, and corners of lateral facets, etc.In order to automatically find the dominant features points in cervical vertebrae, Rochies and Winter proposed researches about detection of anatomical landmarks and dominant points by matching feature sets derived from 2D wavelet and Gabor transform in CT and MRI images [8][9].The proposed methods used the graph matching algorithm to perform a global search, and the similarity of two feature sets was utilized to localize dominant points.However, the method was less adaptive to morphological deformation and the performance on accuracy was not satisfied.Besides, automatic recognition of spine shape is not only used for surgery applications mentioned above, but also an issue of importance of computer-aided diagnosis (CAD).As the growth of the volume of medical images with the progress of medical imaging techniques, it becomes exhausting to inspect all the data in detail manually.Therefore, CAD system is introduced to improve medical automation and the medical data can be pre-processed automatically and then provide information for assisting diagnosis.The anatomical structure around the spine such as soft tissues, muscles, and glands can be regarded as a planar adjacent anatomical space, so the relative location of anatomical structure is able to be inferred according to anatomical knowledge [10].
The anatomical landmarks detected in images are able to be set as the starting points of image segmentation for automatic diagnosis and navigation.In addition, the points of dominant features in 2D slices can also be seeds for shape modeling, 3D reconstruction and registration.In neck CT images, the cervical vertebrae are significant landmarks for medical application but to automatically extract the precise feature points from the complex images is still a challenging task.In this paper, we propose a method to automatically find the feature points of cervical spines in CT slices as shown in Fig 1 based on geometric analysis in companion with anatomical knowledge.
In Section 2, the geodesic active contour model referencing the gradient information is mentioned to extract the cervical spine.We moreover propose a method to automatically set the initial conditions including initial contour and initial distance field for active contour iteration.Then in Section 3, the shape contour extracted from the previous step is coarsely segmented and modeled with a proper scale by means of the curvature scale space (CSS) analysis, and multi-scale geometric analysis is then applied to identify the dominant feature points at each segment.As shown in Fig. 1, the dominant points proposed to extract include the vertices at both sides of vertebral body, near transverse foramens and pedicles (points 1 and 2), the corners of the facets (points 3 and 4) and the corners of spinous process.Furthermore, the vertebral body and lamina regions can be further inferred based on the four determined dominant points.In Section 4, the experiment is carried out on 250 neck CT slices, and the points found out are finally examined with www.ijacsa.thesai.org the points pointed out by the clinical experts to evaluate the success rate and accuracy.In Section 5, we discuss the result and conclude the work.
Typical cervical spine shape with dominant features labeled.

II. VERTEBRAE EXTRACTION
As shown in Fig. 2, the cervical spine locates at the center of the neck and appears in high brightness, and there are various textured soft tissues around the cervical spine.The air path in relatively low brightness is adjacent to the cervical spine.
A cervical CT image and the relative location of the air path and the cervical spine.
In order to extract the spines in CT images, there have been several segmentation methods proposed, including modelbased segmentation, adaptive thresholding, multi-scale canny edge detection and active contour algorithm [1,[10][11][12][13][14][15].Among these methods, gradient-based methods are considered to perform better accuracy than gray level thresholding because the magnitude and direction of the gradient can be used to accurately locate the edges.Besides, in order to segment deformable objects, active contour methods are considered as an effective method to generate continuous boundaries.Geodesic active contour (GAC) is an active contour model (ACM) based on the relation between active contour and the computation of geodesic or minimum distance field [16][17].The initial contour deforms and gradually converges toward the region boundaries through iterations controlled by the gradientbased stopping function with updating the distance field [16].In many applications of medical image segmentation, the initial contours were manually placed near the targets for more effective converging properties and less computation.However, manual placement is not appropriate for automatic segmentation.In this section, we describe a method about automatic placement of the initial contours for delineating the cervical spines and the initial distance field for GAC computation.The method is summarily shown as the schematic diagram in Fig 3.

A. Geodesic Active Contour
The main idea of the active contour model is to minimize the energy associated to the contour derived from the image gradient and the contour curvature.In order to formulate the energy for GAC computation, a contour C in an image is represented by the parametric vector equation: Hence, the energy function of the GAC contour model comprising the internal and external energy terms can be described as the terms below: t is the arc length parameter, E int is the internal energy while E ext is the external energy of contour C. Let C 0 be the initial contour for active contour iterations and g(x) denotes a monotonically stopping function which conducts the contour converge toward the boundary points based on the direction and magnitude of gradient.In order to deform an initial contour www.ijacsa.thesai.orgtowards local minima points of the energy function in the image, the steady state solution is given by where κ represents the curvature derived from the equation: where ．means the derivative, i.e., dt and N  is the unit inward normal.The curve evolves according to the steepest-descent method to deform the initial curve C 0 based on the curvature and gradient.This geodesic problem can be solved by introducing the level-sets approach [16][17] .In level set formulation, the contour C is regarded as the zero level-set of a function, so a contour can be represented as a distance map measuring the minimum distance from a point to the contour.Therefore, the curve evolution ( 4) can be represented by where u is a signed distance field of a contour C, and C can be regarded as the zero level-set in u with u 0 denoting the initial distance field.The level set method evolves a contour by propagating the wave front.The fronts move ahead with a velocity V and arrival time T, and the level set front propagation equation is given by [18] The distance map u is iteratively updated by means of computing the narrow band near the existing front and solving the propagation equation to bring new pixels into the narrow band.The curve evolution operation keeps until the front does not move or the number of iterations approaches the limit.However, in practice, GAC iterations require an initial contour C 0 in companion the corresponding distance field u 0 , and the initial condition significantly affects the result.

B. Automatic cervical vertebrae extraction
In order to obtain C 0 and u 0 for GAC, we start with thresholding the original CT images and the region of gray level higher than the threshold is set to 1 or 0 otherwise.Because the spine regions in a CT image appear in uniformly higher brightness than other regions, the threshold should be chosen higher than the result from Otsu thresholding method [19].However, other regions with high brightness as the spine region would possibly be found, like mandible bones and carotids.
Fortunately, the air path shows distinguishable darkness in cervical CT images and anatomically locates nearby the cervical spine.Therefore, after the morphological operation is proceeded to remove the regions of small area in the binary image, the large blob nearest to the air path is selected.Erosion operation with a 3 3 mask is then involved in to extract the external boundary by collecting the removed elements.The approximate contour of a spine as a result could be sketched out and arranged into a series.The initial contour can be set closed to the real boundary, so it is an effective initial contour for active contour iteration.Different from traditional methods which simply place a circle or a square around the target, the initial distance field can be generated based on geometric relation.In order to build the distance field with an arbitrary contour C 0 , the field need to be derived by distance transform which measured the minimum distance from a point to C 0 .The distance transform can be determined by ) , (

III. SHAPE ANALYSIS AND FEATURE IDETIFICATION
In this section, we describe a coarse-to-fine method to deal with the deformable shape.The spine shape contour is at first coarsely segmented and modeled with CSS, and the dominant feature points on the segmented curves are then figured out by analyzing the detail features in fine scales.The method in this section is demonstrated in the summarizing diagram in Fig. 5.

A. Coarse contour segmentation with CSS
After shape contours are extracted from the original images, the dominant points on the contour which are meaningful for deformable shape recognition need to be figured out.For convenience and symmetry, the point, which is closest to the center of air path, is assigned as the starting point of the closed contour.For shape registration and recognition, geometric points invariant over rotation, scaling and partial occlusion are in considerable importance for the deformable shapes.In general, curvature is a significant property of curves, and the local maximum points or zero-crossing points of curvature profile are considered as meaningful points of the shape [20][21][22].However, local maximum points of the curvature, which could correspond to the corners or vertices of the contour, are too sensitive and easily affected by noise.Because zerocrossing points demonstrate the intersection between a concave contour segment and the adjacent convex segment, as the scale gets higher, the neighboring zero-crossing points also gradually merge together.
Zero-crossing points of the curvature are more adaptive geometric features for deformable shape analysis.CSS is a multi-scale method of collecting zero-crossing points of curvature of a closed contour derived from each scale has been proven as an effective method for shape description and matching over scaling, rotation, partial occlusion and deformation [23][24][25][26].The curvature scale space image (CSSI) is a binary two-dimensional image that records the position of inflection points of the curve convoluted by different-scaled Gaussian filters.In CSSI, along the horizontal axis is the normalized arc length of the contour from 0 to 1; along the vertical axis is the scale parameter.As the standard deviation of Gaussian functions varies from small to large, the contour is gradually blurred while details are gradually eliminated.The multi-scale curvature can be computed by the following equations. where where As a result, a contour can be represented by a CSSI image with several CSS contours corresponding to segments of the shape contour in it, and the CSS contours are constructed by zero-crossing points at different scales.Each CSS contour in the CSSI represents a concave or convex segment of the corresponding shape contour.In the CSSI derived from a typical cervical spine shape, four significant contours of the largest σ are depicted and correspond to four segments of a spine shape contour.With the other four segments squeezed between every two of the fours corresponding to the four significant contours in CSSI, the spine shape contour can be mainly segmented into eight periods of curve.Because a zerocrossing point is regarded as a breakpoint of a concave segment and another convex one, a shape contour can be divided into several meaningful segments by localizing the zero-crossing points.
Each shape contour has its peculiar arrangement of zerocrossing points.Ming et al. proposed a CSS-based method for pattern matching by comparing the CSSI line by line [23], because the zero-crossing points at each scale can be recognizable features of shape contours.In order to effectively separate the shape contour of spine into the eight main segments, an appropriate standard deviation value of Gaussian filter σ, which is related to scale needs to be chosen.Let σ n denote the peak of the CSS contour having the nth highest σ in CSSI.The threshold for choosing the analysis scale in this paper is determined by  Fig. 10 demonstrates a spine shape model based on CSS.It can be observed that a spine shape contour can be segmented into several CSS contours, and there are four main CSS contours labeled as cA 1 , cA 2 , cA 2 and cB 2 in two symmetric pairs in a CSSI of a spine shape contour.cA 1 , cA 2 , cA 2 and cB 2 are corresponding to the four contour segments A 1 , A 2 , A 2 and B 2 as shown in Fig. 10(b) respectively.The segment corresponding to the spinous process locates at the period labeled as C B between the pair of CSS contours (cB 1 , cB 2 ) near the middle of the horizontal axis.Also, the segment C A corresponding to the vertebral body can be deduced by referencing another two apparent CSS contours (cA 1 , cA 2 ) at the both sides of the starting (end) point.The segments squeezed by (cA 1 , cB 1 ) and (cA 2 , cB 2 ), are corresponding to the facets at both lateral sides of the spine, are labeled as C 1 and C 2 respectively.Eight zero-crossing points among (cA 1 , cA 2 ) and (cB 1 , cB 2 ) can be extracted at scale σ analysis , and the eight points are used to separate the spine contour into eight segments.The apparent segments corresponding to the vertebral body, the lateral facets and the spinous process, as a result can be coarsely indicated in CSSI.
In order to find the two main symmetric pairs from the original CSSI as Fig. 6(e), the symmetry property is utilized.(cA 1 , cB 1 ) and (cA 2 , cB 2 ) are symmetric against the starting point of contour, which is contour point closest to the air path center at C A .From the CSS contours with the highest σ, every four CSS contours such as (cA 1 , cA 2 ) and (cB 1 , cB 2 ) in Fig. 10 are extracted at a time.The contour segments within each two CSS contour pairs are extracted, such as the segments corresponding to the periods of (cA 1 , C 1 , cB 1 ) and (cB 2 , C 2 , cA 2 ).The extracted segments are measured for their curve similarity, and high similarity means high symmetry of the contour segments.The four CSS contours in CSSI corresponding to the two contour segments having the highest curve similarity are considered as the four main CSS contours of a spine shape and labeled based on the model.The similarity of two curves is estimated by measuring the difference the two curves.Let f 1 and f 2 be two curves with N uniform sampling points

 
)) ( ), ( ( where t 1 and t 2 are integers, (x 1 , y 1 ) and (x 2 , y 2 ) are the points belonging to f 1 and f 2 and . The dissimilarity D of two curves is evaluated by measuring the distance of the points between two curves.

B. Shape analysis and dominant point identification
After coarse segmentation, dominant points or landmarks can be finely determined by analyzing the contour based on geometrical properties of each segment.Curvature is an important property of a curve, and has been widely used for precisely finding vertices of shape.However, it is so sensitive to noise and small variation of the contour that it is not appropriate for deformable objects or objects in complex images like cervical CT slices.Turning angle is another useful geometric property of curves for comprehending the local variation [27][28][29].Xu et al. applied the turning angle property for curve evolution in automatic spine shape analysis [30].The bending angle of two adjacent line segments is computed and normalized for evaluating the contribution to the whole shape."Included angle" was defined in [1] to figure out the sharp convex characteristics between two adjacent elements.In our research, to evaluate the curve segments, bending angle profile along the whole contour is calculated by the following equation: www.ijacsa.thesai.org ) ) , ( where where  means the inner product operator,  means the norm of a vector and d is the scale parameter.The value of turning angle is closed to but not larger 180 in degrees at the period without acute variation.On the contrary, vertices may locate at the points of local extreme value.As the scale parameter d grows from small to large, only salient vertices can be preserved and more noise is eliminated.Therefore, a multiscale method is introduced in this work.The turning angle profile of a contour is at first derived at a higher scale d 1 , and the local minimum points recognized as dominant points are coarsely localized.Then, the points are adjusted by referencing the nearest minimum points derived at a finer scale d 2 , as indicated in Fig. 12.For suppressing undesired noise, the turning angle profiles convolve with a Gaussian function before extracting the local minimum points.However, the vertices can not be directly discriminated as a concave or a convex one by  but it is necessary for verifying the dominant shape features at each segment.After points of local extreme are localized, each extracted point is recomputed for it curvature by (11) at proper scale associated with d 1 .If the curvature is positive, it means that the angle is convex, or concave otherwise.As shown in Fig. 1, points 1 and 2 are two concave vertices and the segment within point 1 and point 2 indicates the vertebral body.In order to locate the vertices with precision, three adjacent segments of the contour corresponding to the period A 1 , C A and A 2 in Fig. 7 are introduced to compute the geometric information.The two concave vertices result in two local minimum points at both ends of the corresponding segment of turning angle profile.Moreover, the arc segment of the vertebral body can also be accurately determined within these two vertices.Besides, the spinous process is considered to locate at the period C B , which is the middle segment of the whole shape contour.The apparent angles within this period are important because the angles may denote the corners of spinous processes or bifurcation which are apparent landmarks of cervical anatomy.The corners of spinous processes are convex angles, and if bifurcation exists, there will be another concave angle between the two convex corners.Fig. 12 shows the convex and concave vertices in the cervical spine denoting with different marks.In Fig. 1, the facet corners at points 3 and 4 are convex vertices with angle larger than 90 in degree.The "corners" are not only important for nerve root injection operations [2,4], but also for determining the position of lamina periods.In addition to the points labeled in Fig. 1, at the middle of lateral facets C 1 and C 2 , there are sometimes two concave vertices near the facet corners as the points denoted with crosses in Fig. 12.
If existing, the concave vertices are also dominant points near the foramens and spine pedicles as points 1 and 2. The region inside the spine contour between the point 1 or 2 with the adjacent lateral facet concave vertex could be inferred as the spine pedicles denoted with dotted line in Fig. 12.We evaluated the proposed overall framework with two criteria: success rate and accuracy.The success rate was defined as the relative number of dominant points that localized at acceptable positions, and the accuracy was defined by measuring the distance of the points derived from the algorithm with the dominant points drawn by the clinical experts.The accuracy Acc is calculated as follows: where M is the number of points recognized as success.S i is the pixel of the ith point recognized as success derived from our algorithm and R i is the closest reference point drawn by clinical experts.
The average success rate was 93.4%, and the average success rate per vertebra was within the range 70%-100%.The coarse scale of turning angle d 1 could affect the average success rate from 81.7% to 93.4% in our experiment.The reason is that corners with less curvature are blurred at higher scales such that the detailed information is suppressed.The average accuracy is 0.37 mm while the fine scale of turning angle d 2 is 15.

V. DISCUSSION
The experiment results show that the success rate and accuracy are affected by the scale chosen for extract the dominant vertices.The scale relies on the resolution of the CT images, and if the resolution gets higher, the analyzing scale can also get larger.The prominent segments of vertebrae could be adopted as the landmarks for modeling, registration and clinical evaluation.The method we proposed for automatically extract the spine contour can effectively sketch the main contour of the vertebrae without manually setting the initial contour.Unsatisfied results are mostly caused by wrongly contouring, and might occur at the bright blobs connected with the spine without obvious edges.Fortunately, scale-based segmentation could overcome partial deformation.Furthermore, in continuous CT slices, the derived contour can also be the initial contour of the adjacent slices and the distance field for iteration can also be deduced.This is also an extended application of the proposed method in this paper to 3D scene.Turning angle is the main idea for geometric analysis in our study, and dominant points can be recognized by finding the local extreme points of turning angle profile.The vertices with low curvature, such as the facet corners and vertebral body vertices have lower accuracy and success rate than the vertices with sharper angle like spinous processes.The corresponding results derived in adjacent slices could be collected for adjusting, so the performance can be improved.For the purpose of modeling, the dominant points can be applied for building the model mentioned in [14] and [15], and the polygonal approximation can also be deduced.Points 1 and 2 can be used for segmenting the vertebral body, and points 3 and 4 can be used for determine the facet corners.Besides, the dominant points on spinous processes and bifurcation are important landmarks on C7 spine, and can be accurately figured out in this paper.In [1], the accuracy defined by MDCP was 0.14mm with pixel size in-plane ranging from 0.233 to 0.309mm.Comparing with points 3 and 4, which are also landmarks for screw insertion operation, our algorithm performed MDCP accuracy in 0.35mm with pixel size in-plane of 0.78mm.We believe that the detected dominant points are capable for operation assist and the accuracy can also be improved if the image resolution can be higher.

VI. CONCLUSION AND FUTURE WORKS
In this paper, we propose a method for automatically extracting the shape features of cervical vertebrae in CT images.With shape analysis, dominant points of the extracted contour can be figured out.The major contribution of the work is that the proposed method can automatically segment the dominant feature points of shape or landmarks used for operation guiding, therapy planning and model registering.Many proposed models of vertebrae can also be built with more precision by implemented the proposed framework.The framework can also be applied in other the thoracic vertebrae and lumbar vertebrae with adjusting the shape model based on anatomical knowledge.Future works include not only extending the proposed method to other spines, but also building an interactive system for aiding surgery and treatment planning.

Fig. 3 .
Fig. 3. Schematic flow diagram of the proposed method

Fig. 4 .
Fig. 4. (a) An extracted initial contour (b) the corresponding distance field of (a) representing with gray level.The brighter means the point has longer distance to the contour.

Fig. 5 .
Fig.5.The schematic diagram of our algorithm of extracting the dominant feature points by analyzing the shape contour.


is the convolution operator.Function h(t, σ) denotes a zero-mean Gaussian function with kernel size parameter σ and σ is also referred as the scale parameter.If the curve with smoothing parameter has a zero-crossing point at location s on scale σ, we set 1

Fig. 6 .
Fig.6.(a) is the original image.(b) is the extracted contour.(c) is the contour after convolving with the Gaussian filter of σ =15.(d)contour after convolving Gaussian filter of σ =30.(e) the CSSI of the contour in (b)

Fig. 10 .
Fig.10.(a) CSSI of a cervical contour (b) The corresponding segments of contour labeled in (a).

Fig. 11
Fig.11 Measurement of the dissimilarity of two curves.

Fig. 12 .
Fig.12.Turning angle profiles of a cervical spine contour with different scales and the some corresponding vertices at the two scales are pointed out.

Fig. 13 .
Fig.13.(a)(c) Results of automatically contouring.(b)(d).Results of figuring out the dominant points of (a) and (c).○ denotes the dominant points derived from our algorithm and × denotes the dominant points pointed out by the experts.