Automated Determination of a Reference Point for the Diagnosis of Spinal Instability
Introduction
Spinal Instability is a back ailment characterized by an unusual lack of spinal stability and chronic back pain that is caused by the degradation of spinal cord support tissues. Patients who suffer from this ailment require spinal fusion surgery in which two neighboring vertebrae are fused together, thus providing additional spinal stability. The current procedures that are used to determine which vertebrae need to be fused are both difficult and painful. Haughton, et al. [1] has developed a new methodology that will make the same determination simpler and non-invasive. The new procedure relies on a combination of hardware and software that analyzes a set of computer tomography (CT) scans to determine where the vertebral fusion will be most effective to correct Spinal Instability. The software, however, requires a very steep learning curve to make the correct decision. The aim of this project was to design and implement a Matlab-based software program that will drastically reduce the software’s learning curve by automatically selecting necessary points from the CT images.
Before further discussion, it is imperative to understand the physiological background of Spinal Instability and its associated diagnostic techniques. Spinal Instability occurs when the sheath of longitudinal ligaments surrounding the spinal column (Figure 1) begins to degrade. When these ligaments are weakened they lose their capacity to properly hold the vertebrae and intervertebral disks in place. Normal intervertebral disk motion is in the vertical plane with rotation about the horizontal plane (a bending motion). When a person suffers from Spinal Instability, these disks can rotate about the vertical plane as well, causing severe pain.
Figure 1: The anterior anatomical view of the spinal column displays the longitudinal ligaments that provide support for the vertebrae and intervertebral disks – degradation of these ligaments leads to Spinal Instability.
No direct means for determination of the extent of ligament deterioration currently exists; therefore secondary indicators are used to diagnosis Spinal Instability. One of the most common diagnosis methods instructs the patient to stand for a radiograph (X-ray) in a normal upright position. A second radiograph is then taken with the patient bent turning to the right and a third with the patient turning to the left. The physician then uses anatomical markers and analytical geometry to determine how much each disk is moving during side-to-side rotation. These calculations are difficult and less precise due to the loss of spatial information when three-dimensional structures are converted to two-dimensional radiographs. If one or more disks are found to rotate excessively, Spinal Instability may be diagnosed and vertebral fusion is performed. Haughton’s method provides more precise measurements than the procedure described above.
Haughton has proposed the use CT scans to measure the rotation of the intervertebral disks. Through the use of a device that rotates a patient’s lower body 12 degrees in either direction while in the bore of the CT scanner, Haughton can more accurately measure disk rotation. Images are captured while the patient is rotated in each direction and these images are then analyzed using a software package that utilizes digital subtraction techniques to determine the degree of rotation in each vertebra, and hence each intervertebral disk. Two image sets are created, one with the patient rotated to the left, another with the patient rotated to the right. These two images sets are then analyzed using a suite of software packages.
During software analysis, the images are viewed in an external medical imaging program to determine the position of several anatomical landmarks using a Voxel coordinate system. The user is required to find a subset of images that represent the desired vertebrae and their pairs in the opposite image set. Using these images, a reference point on the midline of the anterior wall of the spinal column is manually selected using a special software package that allows radiologists to analyze CT and MR images. The user scans through the image set and manually finds the reference point and records its coordinates. Figure 2 shows an example of the targeted reference point based on expert knowledge. The reference point is then passed to a software program that further analyzes the CT images to determine the extent of vertebral movement between image sets. In a healthy subject, the total amount of vertebral rotation is less than one degree, whereas a subject with Spinal Instability can have over 2.5 degrees of rotation.
The aim of this project was to automate the reference point selection process, as this is a critical part of the diagnosis procedure. If the reference point is not accurately and consistently selected, the digital subtraction software will not register the images properly and therefore will not accurately calculate vertebral rotation. An automated reference point selection procedure, as proposed in this paper, could eliminate a major potential source of error in the Haughton Spinal Instability diagnosis process.
Approach
A highly organized, well-designed plan was required in order to accurately and consistently determine the reference point on the CT images. The images are initially provided in the Digital Imaging and Communications in Medicine (DICOM) format. This format, however, provides a major drawback in that it is difficult to work with in Matlab. After some research, it was determined that the Unix operating system provides a simple image type conversion utility that is suitable for the detection algorithm proposed in this paper. The JPEG image format is suitable for the detection process, thus the following Unix command can be used to convert a DICOM image (foo.dcm) to a JPEG image (foo.jpg):
convert foo.dcm jpg:foo.jpg
An example of a DICOM image can be seen in Figure 2 alongside the same image in JPEG format.
Figure 2: DICOM Image of Lumbar Vertebra and corresponding JPEG conversion (red dot on DICOM image is target point).
Following image format conversion, it is possible to easily manipulate the image and the reference point detection algorithm can begin. The first step in the algorithm is edge detection. Using Matlab’s built-in Sobel edge detection function, it is possible to obtain well-defined images with minimal noise. Through experimentation, it was determined that a Sobel threshold of 0.15 was the most effective in removing noise while providing well-defined edges. An example of the result of the edge detection procedure is found in Figure 3. Of particular note, is the well-defined spinal canal (an extremely important feature in determination of the reference point). To simplify the detection process, a region-of-interest (ROI) is used. The detection algorithm requires that the user-selected ROI includes the spinal canal and the closely surrounding area. An example of a properly selected ROI can be seen in Figure 4. Following selection of the region-of-interest, the algorithm is able to perform more intense processing.
At this stage, the processed image had a great deal of noise in it due to all of the edges present around the spinal canal; therefore, a method for minimizing the number of edges surrounding the spinal canal was required. The first step in this process was to smooth the contours surrounding the spinal canal, thus morphological closing was performed. The closing operation is the product of dilation and then erosion of the image – it isolates the spinal canal as one large open area surrounded by closed noise. Through experimentation, it was determined that a square structure element of size five for dilation and size three for erosion were most effective for this algorithm. The result of dilation and erosion of a typical CT image are shown in Figure 5.
With a now well-defined spinal canal surrounded by a small amount of noise, it is possible to do further computations with little interruption. Once again, edge detection is used to differentiate between the spinal canal and surrounding area. At its current stage, the image consists of only black and white pixels (having values of zero and one); therefore, it is possible to use Matlab’s bwperim function. This function finds the edges defined by changes from one to zero or zero to one. Eight-connectivity was used to ensure that slight imperfections in the spinal canal edge did not cause a break in the detected edges. Figure 6 displays the result of this step in the detection algorithm – there are now three major edges. Of the major edges, the largest is the border of the image and can be disregarded. The second largest edge is the exterior of the region-of-interest selection and is insignificant in terms of further processing. The next largest edge is the spinal canal, which will be used in further processing. Now, the algorithm must extract the image of the spinal canal and determine the appropriate reference point.
Figure 3: Image After Dilation and Erosion
Figure 4: Image After bwperim
Matlab includes the bwlabeln function, which can be used to isolate all individually connected objects and label them with individual symbols. This function assigns each pixel an object symbol, such as five, instead of a pixel intensity value, thus making future processing substantially simpler. For visualization purposes, Matlab’s label2rgb function was implored so that each object was represented by a different color when plotted. The result of object detection can be seen in Figure 7.
Figure 5: RGB Color Mapping of Labels
Once labeled, each object in the image now has an associated set of properties, including object area. Objects in the image are sorted by area, so that the spinal canal can be identified. As explained earlier, the spinal canal object will have the third largest area among all objects in the image. It is now possible to isolate the spinal canal based on its label number. A subimage of the spinal cord edge is created via analysis of the image shown in Figure 7 – if a given pixel contains the correct object symbol, it is assigned a value of one; otherwise the pixel is assigned a value of zero. The result of the subimage process is show in Figure 8, this image contains only the extracted spinal canal image.
Figure 6: Edge of extracted Spinal Canal
As the spinal canal has now been extracted, the only remaining task in the algorithm is to determine the reference point based on the canal image. The spinal canal edge can be approximated by a triangle, of which it is possible to determine three vertices. An estimation of the reference point is then found by calculating the midpoint of the two upper vertices of the triangle. The midpoint provides only an estimation as the posterior edge of the spinal canal is typically not a perfectly straight line (as it is in a triangle). This requires the algorithm to perform an additional calculation to determine the actual position of the reference point. Using the estimation and the third vertex of the triangle, the actual position of the reference point can be found. A line can be drawn from the third vertex through the estimation – the intersection of this line and detected spinal cord edge is the desired reference point. The coordinates of this point is the output of the detection algorithm. Figure 9 displays the original JPEG image as well as the location of the reference point, which is denoted by a superimposed white circle.
Figure 7: Original image with detected reference point (indicated by white dot).
Work Performed
The algorithm discussed above was implemented in Matlab, but a UNIX machine was necessary to perform the conversions of the DICOM images to JPEG images. Also, to facilitate ease of use, a Graphical User Interface (GUI) was developed on top of the algorithm implementation. The GUI allows the user to specify four algorithm parameters: file, Sobel coefficient, dilation structure element size, and erosion structure element size. A screen shot of the GUI can be seen in Figure 10.
To begin processing, the user selects the JPEG version of the CT image the he or she wishes to process using a standard Windows open dialog box. After an image has been selected, the GUI will display both a large size version of the image as well as a thumbnail for future reference. The algorithm parameters can now be set as desired, however, the default values have previously been determined as good values for this algorithm. To begin algorithm execution, the user must select the “Run” button. The large size version of the image is then converted to the edge detected image (as in Figure 3) and the user is able to select the region-of-interest within the image window. After the ROI has been selected, the algorithm will continue execution until the reference point is determined. At this time, the GUI will display the original image with a white superimposed circle to mark the location of the reference point (as in Figure 9). Text boxes just above the “Run” button will contain the coordinates of the reference point. If the user wishes, he or she can now use the “<” and “>” buttons to move through each step of the detection process – a total of eight images can be displayed, one for each major stage in the detection process.
Figure 8: Screen Capture of Graphical User Interface.
Results
The detection algorithm has been tested on seven images of different lumbar vertebra that had all been previously expert classified. The output of each execution was a reference point that was then compared with the expert-chosen reference point for the given CT image. The average Euclidean difference between the expected and actual reference point location was 6.81 pixels, well within an acceptable range. Table 1 displays the results of each trial.
Table1: Expert classified (expected) and algorithm detected (actual) of reference pixel location as well as the Euclidean distance between the actual and expected pixel.
During testing, it was realized that some images need slight alterations to the erosion and dilation parameters to ensure proper detection, but the defaults usually work. In addition, some limitations were realized in terms of the region-of-interest. The user must be sure that at a minimum, the entire spinal canal is selected while also limiting their selection as much as possible. Ideally, the user will select a small polygon surrounding the spinal canal and not select the rest of the vertebra. As the area surrounding the spinal canal in the region-of-interest increases, so too does the potential for error in the algorithm. Beyond these two lessons, the algorithm is consistent and accurate in its detection of the reference point.
Currently, the algorithm is being tested by an associate of Dr. Haughton at the University of Wisconsin Department of Medical Physics. Thus far, the algorithm has received strong reviews.
Discussion
The current diagnosis techniques for Spinal Instability are both complicated for the physician and uncomfortable for the patient. The new methodology proposed by Haughton et al. have given the medical community the opportunity to take a large step forward in terms of diagnosis. This new technique if full of promise, however, several barriers must still be overcome to ensure proper diagnosis. Haughton also proposes that a radiological technician should perform the vast majority of the analysis, while a neuroradiologist finalizes the diagnosis. This project offers a means to simplify a complicated portion of the procedure that also previously required a great deal of diagnostic knowledge.
Although there are a number of further improvements and additions that could be made to the algorithm proposed in this paper, it has been shown that automated detection of the reference point is possible. In the future, the Haughton research group could implement a program that is capable of selecting the best CT images for reference point detection, however, at this stage it was assumed that the radiological technician is able to select the proper images prior to inputting them into the detection algorithm. The algorithm proposed here accurately and consistently detects the appropriate reference point to be passed to the rotation detection software. The proposed methodology is a large step in making the Haughton’s diagnostic technique more feasible for widespread use.
References
- Haughton, V.M. Rogers, B. Meyerand, M.E. Resnick, D.K., “Measuring the axial rotation of lumbar vertebrae in vivo with MR imaging.”, Ajnr: American Journal of Neuroradiology. 23(7):1110-6, August 2002.
- Rogers, B.P. Haughton, V.M. Arfanakis, K. Meyerand, M.E., “Application of image registration to measurement of intervertebral rotation in the lumbar spine.” Magnetic Resonance in Medicine. 48(6):1072-5, December 2002.
Page 1 of 10