Segmentation of anatomical structures on chest radiographs

1

Segmentation of anatomical structures on chest radiographs S. Juhász(1) ), Á. Horváth(1), L. Nikházy(1), G. Horváth(1), Á. Horváth(2) (1)

Budapest University of Technology and Economics/Department of Measurement and Information Systems, Budapest, Hungary (2) Innomed Medical Co. Budapest, Hungary

Abstract— In this paper we present a solution for segmenting anatomical structures on chest radiographs. First we show an algorithm for the lung contour detection, then we describe a method for finding the ribs and clavicles. The results of these procedures are used as input for a bone shadow elimination algorithm. Different implementation results are also discussed, like C++ and a parallel solution on GPU. Keywords— radiograph, chest, segmentation, ASM, dynamic programming.

I. INTRODUCTION

Recently intensive research and development work has been started throughout the world to develop Computer Aided Detection or Diagnostic (CAD) systems to help pulmonologists/radiologists in the evaluation of chest radiographs, where CAD systems usually serve as second readers in the evaluation of radiographs. Approximately two years ago a Hungarian consortium has also started a research/development work to extend a recently developed PACS with CAD functionality. The general goals and the first results of the system were presented last year in the World Congress of Medical Physics and Biomedical Engineering, Munich [1]. Computer aided evaluation of X-ray radiographs needs complex image processing/pattern recognition algorithms where first the images should be preprocessed, and the detection of abnormalities is done in the preprocessed images. Preprocessing means that the contours of the lung fields and the heart, as well as the contours of the bones – the clavicles and the rib cage – should be determined. The findings of these contours may serve two goals: (1) the shape of these anatomical parts, especially the shape of the lung fields may have diagnostic meaning, (2) having determined the contours of the bones and the heart, there is a chance to eliminate the shadows of these parts, “cleaning” the whole area of the lung fields from the “anatomical noise”, and making possible to “look behind” these parts. The suppression of the shadows of “disturbing anatomical parts” may significantly improve the performance of nodule detection, and may help in reducing false positive hits. The goal of this article is to present several preprocessing steps: lung contour detection, rib and clavicle contour detec-

ARON_Sanyi_ver2_corr.doc

tion, bone shadow elimination. Then we briefly discuss the results of an efficient GPU-based implementation.

II.

LUNG CONTOUR DETECTION

A. Introduction A typical first step of chest X-ray image analysis is the detection of the lung contour. If we have the position of the lungs then the detection of other objects such as the heart, ribs and clavicles will be easier as we already know where to look for them. The size and shape of the lung contour has diagnostic value without further processing too, as it can show cardiomegaly and pneumothorax. We tried several techniques for the segmentation: AAM [1], ASM [2], kNN-classification [3], curve fitting, texture analysis and other ad-hoc methods. So far ASM gave the best and most robust results. Other researchers came to the same conclusion [4]. We define the lung fields as the set of pixels for which the X-rays go through the lung, but without the areas of the heart, mediastinum, aorta and the areas under the diaphragm. On Fig. 1 the top left image shows an example lung contour with this definition. B. Active Shape Mode (ASM)l Active Shape Model is a parametric model of a curve or contour where the parameters are determined from the statistics of many sets of points obtained from different contours of similar images using principal component analysis (PCA) In ASM the boundary of the object is determined by n points. From these points a descriptor vector is created as

x = ( x1 , y1 , x 2 , y 2 ,..., x n , y n ) T

(1)

where xi and yi are the two coordinates of the i-th point of the curve. For the principal component analysis we compute the mean shape of the s training vectors

x= and the covariance matrix

1 s ∑ xi s i =1

(2)

2

S=

1 s ∑ (xi − x)(xi − x)T s − 1 i =1

We choose the first t largest eigenvalues

(2)

λi

of the cova-

riance matrix. As only the most significant eigenvalues are taken into consideration the number of parameters in this parametric description is much less than the number of contour points, n. The corresponding eigenvectors are arranged into the matrix Φ = (φ1 φ2 ...φt ) . The model parameters are computed by

b = Φ T ( x − x)

(3)

and an approximation of the shape from the model parameters is obtained as:

x ≈ x + Φb

(4)

For a given contour a proper b vector is looked for. where all components of b are constrained in the interval ± m λi , where m is a properly selected constant. The search of the parametric description of the contour of an object starts from the mean shape, than two alternating steps are applied until some convergence or a certain number of iterations is reached. In the first step we try to move each contour point perpendicular to the contour. Several positions are tested on both sides of the contour and the best position is chosen. To decide the best position an intensity gradient profile model is built at each contour node and image resolution during the training. To find the best new position for the contour point Mahalanobis distance is used. After all the contour points are updated, to fit a model to the new point set the best b parameter is looked for according to (4). The whole process is repeated several times while the spatial resolution of the image is increased. C. Training and results For obtaining the ASM-based contours a publicly available reference image database of 247 radiographs by the Japanese Society of Radiological Technology [5] has been used. The delineation of these pictures was made by the research group of van Ginneken [6]. We divided the image set to two parts. One half was used for training and the other for testing the algorithm. Several implementations were made. The execution time per image was 20s of the Matlab version and less than 6s of the C implementation. We used one core of a 2.3 GHz Core 2 Duo machine for the test runs. In both cases most of the time was spent on the Gaussian blurring.


The execution time of 6s is still not acceptable in realtime applications, in clinical examinations. Therefore we implemented the contour search algorithm on the graphical processing unit (GPU). This implementation needs fundamentally different programming techniques, because the program runs concurrently on many threads. The nVidia CUDA system was used for this implementation. We managed to reduce the execution time to 250ms. The ASM search algorithm itself took 30ms on average, the rest of the time was spent on Gaussian blurring and memory copying. We used an nVidia GeForce 9800 GTX for testing, which contains 128 processor cores. It is important to note that significantly better results cannot be achieved with a more advanced GPU, since the memory copies between the host computer and the GPU take remarkable part of the current execution time. We evaluated the results on a per-pixel level. We were able to get an average sensitivity of 0.956 at the specificity level of 0.984.

III.

BONE DETECTION

A. Algorithm Looking at chest X-ray images an apparent regularity can be observed at first glance: the ribs have a quasi parallel arrangement. So it seems feasible to model the ribcage somehow. Even though the two sides of the ribcage appear to be symmetric they have to be processed separately. The possible gain from utilizing this symmetry isn't worth while considering the complexity of handling the occasional disorders and irregularities. To build a complete model which is capable of enumerating the ribs one-by-one and to specify their positions would require too many parameters and still wouldn’t be accurate enough to handle all the anomalies. Thus we chose a different approach. We restrict our model only to the slopes of the ribs. Our model assigns a slope to every point of the image, this way it neglects the position information. These slopes can be described by a function with two arguments and graphed by a slope field. Several functions were investigated for this purpose and it is revealed that 3 parameters are enough to describe the ribcage. The final version of the function was produced by applying principal component analysis on a third order rational function of two variables. On Fig. 1 the top right image shows an example of this model. During model fitting the missing coefficients are found by an exhaustive search. In order to measure the fitness of the model to the image an edge detection step has to be performed. After selecting some points from the area of the lung, the best fitting arcs could be found in the vicinity of

3

each point, resulting in a series of arcs defined by their positions, slopes, and curvatures. This kind of edge detection assures that the edge segments have a uniform distribution over the lung. The position and slope parameters of the arcs can serve as samples to adjust the model. From the obtained model quasi parallel curves can be generated including the real rib outlines. The lower and upper outlines of the ribs are not discriminated. By selecting the best fitting curves we can get to the approximate outlines of the ribs. This selection process is not straightforward, because next to a good match several other curves are present with similar high fitness values. Hence a local maximum selection has to be incorporated which considers the admissible distances between these outlines too. If these restrictions are well configured we can get the most likely approximation of the outlines by applying a dynamic programming on these constraints. To refine the approximate curves we introduced a parallel version of the dynamic programming active contour algorithm. This method handles the lower and upper outlines of the ribs simultaneously this way ensuring the quasi parallelism of these curves. Not only parallelism but other constraints can be incorporated: the rib thickness also can be restricted to the admissible interval and this method can foster the smoothness as well. By applying these constraints false results can be strongly reduced. This algorithm takes the middle-curve as a basis and then shifts its points along the normal of the curve at each point, this way creating several copies running parallel to the starting curve. Then these copies are split into short segments. The segments can be arranged into a matrix where each curve represents a row and each curve-segment can be assigned to an element in this row. The cell values represent the fitness of a given segment to the image. Now the goal is to find two quasi parallel path between the left and right side of the matrix with maximal sum. From every column two cells have to be selected. A pair of paths can be described by a mid path and a thickness value in every column. Thus the problem is mapped to another domain. The thickness defines the distance between the two selected cells in the given column. Constraints can be introduced for the adjacent path elements and for the adjacent thickness values. The cell which represents the middle of the rib cannot move more than one row at every column border. The adjacent thickness values are also restricted to change only one unit at a time. Constraints for the overall thickness of a rib can also be applied here. This problem can be transformed to a longest path problem in a graph, where nodes represent the cell and thickness pairs and the edges express the allowed steps between the adjacent columns. The edge values reflect the value of the cell which are pointed by the edge plus the additional penal-


ty values for changing row or changing thickness value. By solving this problem the rib outlines can be obtained with decent accuracy as seen on the middle left image of Fig. 1.

Fig. 1 Results at the main steps of the processing of the first image in the JSRT database. Images at the bottom show a magnified part from this image

This method was applied for the clavicles as well. But in contrast to the ribs’ case it is not required to find approximate clavicles in advance. Predefined clavicle templates can be used instead. If the templates have sheer slope and stark convex curvature they can only fit to clavicles when applying this method. The templates are placed on the images and their relative position to the lung outline was estimated by analyzing several radiographs. This segmentation data is then used to remove the bones from the images in order to enhance the structure of the lungs. The elimination is based on creating intensity profiles on vertically differentiated images, and then subtracting

4

them from this image and returning to the original domain by integration. The bottom images of Fig. 1 depict the result. The execution time of the current algorithm is around 2 minutes on a 2.3 GHz Core 2 Duo machine, which is far from acceptable, but the introduction of GPU accelerated subroutines are already begun successfully. The last part of the procedure which refines the outlines and originally run for 30 seconds could be reduced to 1 second, but including the transfer overhead between the video card and the memory it raised to 6 seconds. The overhead can be mitigated by executing longer continuous parts on the GPU.

results were good enough to be used for further processing steps.

B. Results

The presented clavicle outline detection algorithm performed well on digital radiographs, but further enhancement is still possible to handle the fainter contours of the JSRT images. The rib detection algorithm found enough ribs to produce a clean image of the lung after eliminating the shadows, but in some irregular cases it still skips too many ribs, however there are a lot of possible ways to enhance it.

We tested our algorithm on two sets of images. The first was the first 100 images of the JSRT database [5] and the other was from Innomed Medical Co. The former contains images taken by analog devices and is quite homogeneous. The shadows of the bones are faint but the whole scene is fairly clear. The latter set came from digital machines. It is heterogeneous; it has shoots from different devices with different resolutions and different photon energies. The bone shadows are darker, but the details are blurry and distorted by noise. The patients are mostly elderly having different diseases resulting in abnormal bone and tissue structure. Thus the latter set better represents real usage. Table 1 Distribution of clavicles Found Innomed JSRT

97% 81.5%

Accurate 93% 66.5%

IV. CONCLUSIONS

It was shown that the lung contour detection can be solved with decent accuracy while keeping runtime low enough with the use of GPU for real-time evaluation. The

Percent of cases

JSRT 40

30

30

20

20

10

10

0

Found Correct

0

0

1

2

3

4

5

6

7

8

9

10

0

1

2

3

4

5

6

7

8

9

10

Fig. 2 Distribution of ribs per lung

ACKNOWLEDGMENT This work was partly supported by the National Development Agency under contract KMOP-1.1.1-07/1-20080035.

Wrong 3% 18.5%

We were not able to undertake our system an objective examination, because we do not have reliable information about the position of the bones. Thus the classification of results was done by visual observation. The resulted outlines were divided into three categories. If it is an outline of a real bone than it is "found", if the contour follows the bone generally than it is "correct", and if there is no exactly assignable bone for it than it is "wrong". The number of found, correct and wrong ribs was counted per lungs, because the sides were processed separately. Fig 2 shows the distribution of these aggregates. There were 0.28 and 0.24 wrong outlines per lung respectively. In the case of clavicles these categories were aggregated only on the whole set.


Innomed 40

REFERENCES 1. 2. 3. 4.

5.

6.

Cootes T F, Edwards G J, Taylor C J (2001) Active Appearance Models. IEEE TAMI 23(6):681-685 Cootes T F, Cooper D, Taylor C J, Graham J (1995) Active Shape Models – Their Training and Application. CVIU 61(1):38-59 Goldstein M (1972) k-Nearest Neighbor Classification. IEEE TRANS.. INFORM THEORY 18(5):627-630 van Ginneken B, Frangi A F, Staal J J, tel Haar B M, Viergever M A (2002) Active Shape Model Segmentation with Optimal Features. IEEE TRANS. MED IMAGING 21(8):924-933 Shiraishi J, Katsuragawa S, Ikezoe J, Matsumoto T, Kobayashi T, Komatsu K, Matsui M, Fujita H, Kodera Y, Doi K (2000) Development of a digital image database for chest radiographs with and without a lung nodule: receiver operating characteristic analysis of radiologists’ detection of pulmonary nodules. AM J ROENTGENOL 174:71-74 van Ginneken B, Stegman M B, Loog M (2006) Segmentation of Anatomical Structures in Chest Radiographs using Supervised Methods: A Comparative Study on a Public Database. MED IMAGE ANAL 10:19-40 Author: Institute: Street: City: Country: Email:

Áron Horváth Budapest University of Technology and Economics Magyar tudósok körútja 2. Budapest Hungary [email protected]