A Dynamic Programming Approach for Heart Segmentation in Chest Radiographs Aarti Raheja, Jason Knapp, Chunsheng Fang {araheja,jknapp,vfang}@riverainmedical.com Riverain Technologies 3020 South Tech Boulevard Miamisburg, OH 45342-4860 Abstract hybrid voting scheme was used to combine the three Chest radiographs are the most routinely acquired exams, methods. Shape models, such as the ASM, have the which makes their use for diagnosis cost effective. In this drawback that their fitting routine can get caught in local paper we present a dynamic programming approach for optima [5]. This effect can become quite pronounced when automated heart segmentation on posterior-anterior (PA) applied to images that differ significantly from those used chest radiographs. The goal of the proposed algorithm is to to build the model. This point is particularly important in provide an accurate and reproducible method for heart our application as abnormal hearts are precisely what we’re segmentation, which can then be used to detect certain trying to detect. For this reason, we opted for a different cardiac abnormalities. Our method has several advantages approach. over previous methods, and provides superior performance to previously published results. One important use of heart segmentation is the measurement of the cardiothoracic ratio. The CTR is an important measurement that can imply cardiomegaly Introduction (abnormally large heart) [1]. The CTR is defined as the maximum transverse diameter of the heart, divided by the Heart segmentation in chest radiographs is a maximum internal diameter of the thoracic cage [6] challenging task. One major difficulty in segmenting the Research in to methods for automatic CTR extraction has a heart is the low contrast found in the mediastinum and long history [6]. Later in the paper we show how the CTR diaphragmatic regions. These areas are difficult to can be used for assessing the quality of a heart visualize even by radiologists. Other aspects that make the segmentation. Although the CTR can be computed without problem challenging include: the significant variation in segmenting the heart, segmentation is useful as it can help heart size across patients, the presence of disease in the radiologists validate the result. Figure 1 illustrates the idea. lungs, and poor breadth holds by patients (leading to lower We use an algorithm based on dynamic programming contrast on the heart boundary). Despite the challenges, (DP) to segment the heart. DP, an important algorithm in development of an automated method for heart Artificial Intelligence [7], is used in applications such as segmentation could provide significant clinical value [1]. finding the shortest path within a graph. DP decomposes a Several methods have been proposed [1][2][3] for complicated problem into simpler sub problems; and, based segmenting the heart. Nakamori et al [1] discuss a method on Bellman’s “Principle of Optimality”, the optimal to segment the heart by detecting points along the heart solution to the original problem can be obtained by boundary, which are then fitted using a Fourier shape combining the solutions to each sub problem. model. This method was used in [2] to automatically In the proposed algorithm we formulate the DP sub compute the cardiothoracic ratio (CTR) in 400 problem in an innovative way. The cost matrix is generated radiographs. Out of the 400 radiographs, 20% required using image information assigning minimum cost to the manual intervention. It was also shown in [3] that the heart pixels having heart edge characteristics. The cost matrix is boundaries outlined by four experienced radiologists had a generated in the polar domain since the heart shape is high degree of variability, which is an important result mostly circular. By using this method we allow the shape when considering how to assess automatic methods. to vary in regions where enough information is present, but Van Ginneken et al [4] discuss several approaches to force the shape to be circular in regions of uncertainty. heart segmentation: active appearance model (AAM), In the next sections we describe our algorithm based active shape model (ASM) and pixel classification. The on dynamic programming in detail, followed by a individual methods performed comparably well, though presentation of extensive experimental results and a significantly better performance was obtained when a conclusion. method as discussed in [10]. The average of these two locations, as shown in Figure 3, is used as the end row value to define an approximate bounding box around the heart region. The top row of the bounding box, as shown in Figure 3, is selected as the location where the heart and the left lung first meet. The bounding box column locations, as shown in Figure 4, are the locations along each lung mask that are at a maximum distance from the central column. Figure 1: Chest Radiograph with heart outline showing the maximum internal diameter of the thorax (ID) and maximum transverse diameter of the heart that is the sum of maximum right heart width (MR) and maximum left heart width (ML). Materials and Methods We used the 247 chest radiographs from the JSRT database to test the method. The JSRT database is available publicly and consists of screen-film images digitized with a 0.175mm pixel size and 2048×2048 image size [8]. The heart annotations for this dataset [9] are Figure 3: Landmark points computed using curvature available and were used to evaluate our method. information on the lung masks In Figure 2 a flowchart of the method is shown. This bounding box is used to define a center and a radius around the approximate heart region. The center is selected as the midpoint of the bounding box and the radius is selected as half of the distance between the end column locations. Figure 4: The heart region determined using the lung mask Figure 2: Algorithm Flowchart Polar Transform Region of Interest around the Heart The border of the heart is roughly circular. For this reason we apply a polar transform defined in equations (1)-(3) to the approximate heart region. We first obtain the ribcage mask and the segmented lung masks from the chest radiographs. This is done using a ( ) ( ) (1) method developed by Riverain Technologies. The lung masks are then used to detect locations where the air, heart, and diaphragm intersect as shown in Figure 3. These √ (2) locations are computed based on a curvature detection ( ) (3) enforcing a smoother result. The total cumulative cost where ( ) is the image in the Cartesian coordinate matrix is defined as follows: system and ( ) is the image in the polar coordinate system. ( ) ( ) (5) The polar transform is applied to the image as shown in Figure 5(a) using the center and radius as defined in the ( ) previous section. To ensure all of the heart is included, the ( ) ( ) ( ) radius is multiplied by a factor α. In this paper we selected a α value of 1.5. The polar domain image as shown in (6) Figure 5(b) is used to compute a cost matrix for the purpose of dynamic programming. where T represents the transition cost. The value “s” is the offset between pixels when going from one column to the next. The value of this offset is not allowed to be larger than a specified value, “k”, depending upon the desired path smoothness. The value of k for our experiments was set to 3 pixels. Pixels outside the lung mask, or those having cost values above a maximum acceptable threshold, are set to the maximum cost value as shown in Figure 6. This causes a straight line to be the optimum path for these regions Figure 5.(a): Cartesian system image with center and (circular arc in Cartesian domain). radius marked for conversion in the polar domain Figure 6.(a): Original Cost Image Figure 5.(b): Polar co-ordinate image expressed in terms of radial distance. Figure 6.(b): Cost Image with non-air pixels suppressed Dynamic Programming In dynamic programming, the most important part is Figure 6.(c): Cost image with pixels having cost values above constructing the cost matrix. Each pixel in the cost matrix a maximum acceptable threshold value set to the maximum cost is assigned a local cost, where we use low cost values for value pixels that have characteristics typical of the heart boundary. The local cost is defined as a linear combination Once we obtain the optimal DP solution path, the heart of individual cost images: segmentation is obtained by transforming the path to the Cartesian domain, as shown in Figure 7. (4) where is the cost based on the gradient magnitude, (a) is the weight assigned to , is the cost based on a smoothed gray scale image and is the weight assigned to . The gradient is calculated by computing the derivative along each column (derivative in the radial direction). The gray scale cost term is defined by first computing a nominal value for the heart-lung border. This is done for each column within a smoothed image. These nominal values are then used to measure each pixel’s deviation from the expected border value. Each (b) (c) local cost term is scaled to the unit interval prior to Figure 7: (a) The optimum path obtained from dynamic combining. programing solution is converted into (b) Cartesian co-ordinate Given the local cost matrix, the next step is to compute system to obtain the heart segmentation with (c) some post the cumulative cost. The cumulative cost accounts for processing. both the local and transitional costs. The transitional term weights the cost of going from one pixel to the next. The Some morphological post processing is applied to make transitional cost we use increases with pixel distance, thus the heart shape smooth and convex, see Figure 7 for an example. Experiments Figure 10. The reason this can occur is that the source of low overlap is generally from the mediastinal and sub We carried out two experiments to validate the proposed diaphragmatic regions, which do not influence the method. First, the algorithm output is compared to the transverse diameter of the heart. manual outlines to evaluate the accuracy of the heart Figure 11 shows the only case with a low overlap score segmentation. In a second experiment, we compared CTR that was not due to the mediastinal or sub diaphragmatic values extracted from the algorithm against those extracted regions. The difficulty here is the fusion of the left lung from manual outlines. The specific aim of this experiment and colon. This leads to an inaccurate estimate of the left- was to evaluate if a reliable CTR estimate can be obtained lower landmark intersection location, which results in even with a low overlap score. significant under segmentation. Fortunately, such an The overlap score, Ω, between the manually outlined occurrence is rare and is left as an area for future heart boundary and the output of our method is defined in improvement. equation (7). 20 CTR Relative Percentage Difference 18 (7) 16 Example 1 14 where TP is the true positive area, FP is the false positive 12 area, and FN is the false negative area. 10 Figure 8 illustrates a summary of the overlap scores obtained by our method. 8 6 Example 2 Overlap Score Distribution on 247 images 4 30 2 0 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Overlap Score Number of Images 20 Figure 9: Scatterplot comparing the overlap score with the 15 relative CTR measure 10 5 0 0 0.2 0.4 0.6 0.8 1 Overlap Score Figure 8: Histogram of overlap scores on 247 chest X ray images. Most of the overlap scores concentrate around 0.87, indicating the high accuracy of our method. Figure 10: Example 1 with Figure 11: Example 2 the lowest overlap score of having low overlap score, The CTR values are computed by detecting the internal 0.63 but a more accurate CTR diameter (ID) of the thorax and the transverse diameter of value the heart (TD = MR+ML, Figure 1). Some typical output segmentations are presented in Figure 12. As can be seen, our proposed method captures (8) the actual heart contour fairly accurately in most of the cases. The ID value was derived from the ribcage mask. The TD values were computed using the heart mask derived Discussion from the algorithm output and the manual outlines. A relative difference between the CTR values was An average overlap score of 0.867 ± 0.046 was obtained computed using the above TD and ID values. Figure 9 from the 247 JSRT images. We find that our method shows a scatterplot comparing the overlap score with the produces outputs that are close to the human observer, relative CTR measure. From this plot we can deduce that a while comparing favorably to the other methods discussed good CTR estimate can be obtained even with a low in the survey paper [4]. The overlap scores in Table 1 are overlap score. An example of such a case is shown in for the three hybrid methods discussed in [4]. These hybrid methods make use of multiple methods making them References computationally intensive. In addition, these methods are supervised approaches whose outputs might not extend to [1] N. Nakamori, K. Doi, V. Sabeti, and H. MacMohan, more atypical cases. By comparison, our method is far less “Image feature analysis and computer-aided diagnosis in complex and has the advantage of making very few digital radiography: Automated analysis of sizes of heart assumptions about the shape of the heart. and lung in chest images,” Med Phys., vol.17: pp. 342-350, 1990. 63 [2]N. Nakamori, K. Doi, H. MacMohan, Y.Sasaki, and S. 3 Montner, “Effect of heart-size parameters computed from digital chest radiographs on detection of cardiomegaly: Potential usefulness for computer-aided diagnosis,” Investigat. Radiol., vol. 26: pp. 546-550, 1991 [3] R. Kruger, J. Townes, D. Hall, S. Dwyer, S. 65 80 93 Lodwick. “Automated radiographic diagnosis via feature extraction and classification of cardiac size and shape descriptors,” IEEE transaction on Biomedical Engineering., vol. BME-19:pp. 174-186, 1972 [4] B. Van Ginneken, M. Stegmann, M. Loog. “Segmentation of anatomical structures in chest radiographs using supervised methods: a comparative study on a public database”, 2004. 69 78 [5] T. F. Cootes, C. J. Taylor, D. Cooper, and J. Graham. “Active shape models – their training and application,” Computer Vision and Image Understanding., vol. 61(1):pp. 38–59, 1995. [6] H. Becker, W. Nettleton, P. Meyers, J. Sweeney, Jr CM Nice. “Digital computer determination of a medical Figure 12: Heart Segmentation. Blue represents user annotation and red represents output of the current method. diagnostic index directly from chest X-ray images,” IEEE Transaction on Biomedical Engineering., vol. BME-11:pp. 67-72, 1964. Conclusion [7] Stuart Russell, Artificial Intelligence: A Modern We presented an algorithm for segmenting the heart region Approach . using dynamic programming. The proposed algorithm [8] J. Shiraishi, S. Katsuragawa, J. Ikezoe, T. Matsumoto, provided an accurate and reproducible method for heart T. Kobayashi, K. Komats, M. Matsui, H. Fujita, Y. Kodera, segmentation. The presented method makes few and K. Doi. “Development of a digital image database for assumptions about the heart shape, has a simple chest radiographs with and without a lung nodule: Receiver implementation, and provides superior performance to operating characteristic analysis of radiologists' detection previously published results. of pulmonary nodules”. AJR., vol. 174:pp. 71-74, 2000. Future work will involve the collection of more data, [9] Image Sciences Institute Research Databases. which is needed for further evaluation and the http://www.isi.uu.nl/Research/Databases/. development of strategies for handling outlier cases. Also, [10] S. Muhammad, M Asif and M. R. Asin. “A new additional image features for improving the local cost term approach to corner detection,” Computer imaging and will be explored. vision, vol. 32:pp. 528-533, 2006. Table 1: Overlap score results compared to a human observer and various methods discussed in [4]. Heart µ±σ min Q1 median Q3 max Human Observer 0.878±0.054 0.571 0.843 0.888 0.916 0.965 Dynamic Programming 0.867±0.0460 0.636 0.846 0.875 0.898 0.944 Hybrid Voting 0.860±0.056 0.651 0.833 0.870 0.900 0.959 Hybrid ASM/PC 0.836±0.082 0.430 0.804 0.855 0.889 0.948 Hybrid AAM/PC 0.827±0.084 0.499 0.791 0.846 0.888 0.957