Organ Pose Distribution Model and an MAP Framework for Automated Abdominal Multi-Organ Localization Xiaofeng Liu, Marius George Linguraru, Jianhua Yao, Ronald M. Summers Imaging Biomarkers and Computer Aided Diagnosis Laboratory, Radiology and Imaging Sciences, Clinical Center, National Institutes of Health, Bethesda, MD 20892

Abstract. Abdominal organ localization is required as an initialization step for most automated abdominal organ analysis tasks, i.e. segmentation, registration, and computer aided-diagnosis. Automated abdominal organ localization is difficult because of the large variability of organ shapes, similar appearances of different organs in images, and organs in close proximity to each other. Previous methods predicted only the organ locations, but not the full organ poses including additionally sizes and orientations. Thus they were often not accurate enough to initialize other image analysis tasks. In this work we proposed a maximum a posteriori (MAP) framework to estimate the poses of multiple abdominal organs from non-contrast CT images. A novel organ pose distribution model is proposed to model the organ poses and limit the search space. Additionally the method uses probabilistic atlases for organ shapes, and Gaussian mixture models for organ intensity profile. An MAP problem is then formulated and solved for organ poses. The method was applied for the localization of liver, left and right kidneys, spleen, and pancreas, and showed promising results, especially on liver and spleen (with mean location and orientation errors under 5.3 mm and 7 degrees respectively). Keywords: Organ localization, maximum a posteriori, probabilistic atlas, pose distribution model.

1

Introduction

The volumes and shapes of abdominal organs can be indicators of disorders, and computed tomography (CT) is commonly adopted for abdominal diagnosis and preoperative planning and guidance. In computer-aided diagnosis (CAD), the identification and segmentation of abdominal organs are essential for further assessment. As in clinical practice, there is a strong CAD incentive for the automated simultaneous detection and analysis of multiple organs, which benefits from interorgan spatial relationship and interaction, for comprehensive diagnosis. Many methods have been proposed for the segmentation of individual abdominal organs from contrast-enhanced CT images especially for the liver [1,2,3,4], and recently the simultaneous segmentation of multiple organs [5,6,7]. Most of these methods relied on prior knowledge of the organs, for example probabilistic atlases

2

[8,9], which are sensitive to initialization/registration, and active shape models [1,4]. The knowledge was utilized to provide an initial segmentation and then refined by fitting the shape model and/or using geodesic or parametric deformable models [4]. In addition, multi-dimensional contrast-enhanced CT data were also employed in appearance-based segmentation using component analysis in a Bayesian framework [10], or using a 4D convolution constrained by a historic model of abdominal soft tissue enhancement [6]. Abdominal multi-organ segmentation remains a challenging task because the sizes, shapes and locations of the organs vary significantly in different subjects. Moreover, these organs have similar appearance in CT images, especially non-contrast data, and are in close proximity to each other. Thus the successful segmentation requires a good initial identification and localization of individual organs, generally performed interactively [5,7,8]. Correct organ localization can also benefit other image processing tasks, including registration and computer-aided detection. Among the most notable automated localization techniques for abdominal organs, Okada et al. [1] initialized the liver segmentation by estimating the abdominal cavity, but it is not certain how well this approach works for smaller organs, e.g., kidneys and pancreas. Yao and Summers [12] used a statistical location model, but the method was limited to estimating only the organ locations without considering the orientations and sizes. Yao et al. [13] simultaneously detected multi-organ locations by finding bounding boxes using principal component analysis and a probabilistic atlas. Due to the large variability of abdominal organ sizes and orientations, however, the location alone cannot completely localize the organs in the abdomen, and thus is not sufficient to accurately initialize other image analysis tasks. Seifert et al. [11] estimated the organ location, orientation, and size using automatically detected anatomical landmarks, semantics and machine learning techniques, but the technical details of the method are not clear. In this paper, we propose a maximum a posterior (MAP) framework for automated abdominal multi-organ localization. Our method finds the poses of multiple abdominal organs, which include not only the locations, but also the orientations and sizes. The paper first introduces a new abdominal cavity normalization to reduce the variability of organ pose caused by different abdominal sizes across subjects. Next, from a training data set we compute a novel pose distribution model that estimates the probability density functions of organ poses using Parzen windows. The method also uses a probabilistic atlas to model the organ shapes, and a Gaussian mixture model to compute the organ intensity profiles. Last, we formulate an MAP problem and solve it for multi-organ poses from abdominal CT images. The method was applied to 12 data sets (5 organs/set) with promising results.

2

Method

Our method estimates the poses of abdominal organs from non-contrast CT images using a maximum a posteriori (MAP) framework. In this work we focus on five organs, i.e., liver, spleen, left and right kidneys, and pancreas. From a group of N training images {I n } nN=1 , these organs are manually segmented and statistically modeled by building a pose distribution model (OPDM), a probabilistic atlas (PA),

3

and a probabilistic intensity profile (IP). For a given subject image, these abdominal organs are then localized using an MAP framework based on these pre-computed statistical models. The method is detailed as follows. 2.1

The MAP Framework

For a given abdominal subject CT image S , the organs are localized by finding the pose parameters that maximize the a posteriori probability, i.e., Θ ( j )* = arg max p(Θ ( j ) | S, O ( j ) ), ( j) Θ

(1)

where O ( j ) is the jth organ, and Θ ( j ) is the pose of O ( j ) in the given image. The organ pose Θ ( j ) is defined using 9 parameters, which include the location c = [c x , c y , c z ] , orientation v = [v x , v y , v z ] , and scaling s = [ s x , s y , s z ] .Using Bayes’ theorem, the a posteriori probability is re-written as p (Θ ( j ) | S, O ( j ) ) =

p (S | Θ ( j ) , O ( j ) ) p(Θ ( j ) | O ( j ) ) . p (S | O ( j ) )

(2)

The denominator in (3) is not related to Θ ( j ) and thus can be ignored. The modeling of the prior p (Θ( j ) | O ( j ) ) and the conditional probability p (S | Θ( j ) , O ( j ) ) is explained as follows. 2.2

Abdomen Normalization

To reduce the organ pose variances, the abdomens of different subjects are normalized to account for the shape and size differences. For this, the vertebrae and the ribs are automatically segmented and identified from the CT scans using the method in [12]. Here we consider four vertebrae, T11, T12, Fig. 1 Illustration of segmented vertebras and ribs, L1, and L2, because they span over and the abdominal cavity bounding box. the location of the majority of the abdominal organs. An abdominal bounding box defines the abdominal cavity and is determined as such: one edge is parallel to the vertebrae line defined by the centers of the four vertebrae and has a length that covers exactly the four vertebrae, and the other edges are defined by finding the minimum bounding box that contains the ribs. The segmented vertebrae and ribs, and the bounding box are illustrated in Fig 1. For normalization, one standard image (denoted as J 0 ) is chosen from the training images, and its abdominal cavity is considered as the normalized space. All the other images are then normalized to J 0 by aligning the abdominal bounding boxes and are denoted as J n for n = 1,..., N − 1 .

4

2.3

The Prior: Organ Pose Distribution Model

The prior p (Θ( j ) | O ( j ) ) is modeled using Parzen window [14] technique from the training data sets, and we call it organ pose distribution model (OPDM). The OPDM for each abdominal organ describes the statistical variability of the organ pose, and is built independently. The pose of organ O ( j ) in the standard image J 0 is defined as the reference pose, i.e., Θ (0 j ) = [c (0 j ) , v (0 j ) , s (0 j ) ] with c(0j ) being the center of gravity, the organ orientations in the three dimensions being v (0 j ) = [0,0,0] , and scales being s (0 j ) = [1,1,1] . The organ poses in the other training images relative to the reference pose are computed using a 9-parameter linear registration on the manual segmentations. For each training image J n the manually segmented organ is registered to the manual segmentation in J 0 to find the transformation T ( j ) (⋅) with parameters [u, r, t], where u, r, and t are the scaling factors, rotation angles, and translation vector respectively. Different organs are registered separately. The organ pose θ n( j ) in J n is then calculated as:

c (n j ) = Tn( j ) (c (n j ) ) , v (n j ) = r , and s (n j ) = t for n = 1,..., N − 1 . The 9 pose parameters are assumed independently distributed. Though these parameters may be weakly related, this assumption greatly simplifies the model and reduces the dimensionality of the OPDM. The probability density function of each pose parameter is estimated using Parzen window method. Let θ k be any of the 9 pose parameters, and θ k( ,jn) be its value in the nth image for organ O ( j ) , then

p (θ k | O ( j ) ) =

1 N

− 1 e 2π h

N −1 n =0

(θ k −θ k( ,jn) ) 2 2 h2

,

(3)

where h is the bandwidth and is estimated using the standard deviation of the data. The OPDM is then constructed as

p (Θ ( j) | O ( j ) ) =

9

∏

p (θ k | O ( j ) ) .

(4)

k =1

2.4

Probabilistic Atlas (PA) and Intensity Profile (IP)

To model the conditional probability p (S | Θ( j ) , O ( j ) ) , a probabilistic atlas (PA) and an intensity profile (IP) are built for each organ. For organ O ( j ) , the PA is constructed using the linear registration results (see Section 2.3) on the manually segmented training images. The probabilistic value on each voxel x i in the image volume is computed as the number of training images that x i is labeled as part of O ( j ) divided by the total number of training sets N. We denote the PA of O j as p (x i | O ( j ) ) . The organ intensity profile (IP), denoted as p (u | O ( j ) ) , describes the probability that a voxel in the organ O ( j ) takes an intensity value of u. It is constructed from the histogram of the Hounsfield unit (HU) value inside the manually segmented organ,

5

and is fitted with a Gaussian mixture model. In the method we use 3 Gaussian components because it produces sufficiently good histogram fitting with low computation. 2.5

The Conditional Probability and Cross Entropy

For simplification, it is assumed that the HU values of voxels inside organs are independent and identically distributed with a known pdf defined as p (u | O ( j ) ) . Thus, the conditional probability in Eqn. (2) can be written as p (S | Θ ( j ) , O ( j ) ) =

∏ p (x

i

| Θ ( j ) , O ( j ) ) p (u (x i ) | Θ ( j ) , O ( j ) ),

(5)

x i∈V

where V is the set of voxels in the image volume, and u (x i ) is the HU value at x i for given pose Θ ( j ) ( p (u (x i ) | Θ ( j ) , O ( j ) ) is not dependent on Θ( j ) ). Let the range of HU values that a voxel can take be [u1 , u 2 ,..., u M ] . For a given pose Θ( j ) , we define the conditional histogram h(u m | Θ ( j ) , O ( j ) ) as the summation of the probability that a voxel belongs to Θ( j ) over all voxels that take an image value of u m , i.e., h(u m | Θ ( j ) , O ( j ) ) =

f (u m , x i ) p(x i | Θ ( j ) , O ( j ) ), with f (u m , x i ) = xi∈V

1, if u (x i ) = u m 0, otherwise

.

(6)

Thus, Eqn. (5) can be re-written as p (S | Θ ( j ) , O j ) =

M

∏ ( p(u

m

| O ( j ) )) h(u m |Θ

( j)

,O ( j ) )

m =1

= exp

M

h(u m | Θ ( j ) , O ( j ) ) log p (u m | O ( j ) .

m =1

(7)

The exponent in Eqn. (7) is the negative cross entropy between the two probability functions h(um | Θ ( j ) , O ( j ) ) and 2.6

p (um | O ( j ) ) .

Pose Estimation using the MAP

The MAP problem in (1) and (2) can be simplified using Eqns. (4) and (7). As the denominator in Eqn. (2) is not related to Θ ( j ) , it can be ignored from the MAP estimation. In addition, the maximization of Eqn. (1) is equivalent to maximizing the logarithm of the a posteriori probability (see Eqn. (2)). Thus the MAP problem in Eqn. (1) is equivalent to solving Θ ( j )* = arg max − H (h(u m | Θ ( j ) , O ( j ) ), p (u m | Θ ( j ) , O ( j ) )) + ( j) Θ

9 k =1

log p(θ k | O ( j ) ) ,

(8)

with H (⋅ , ⋅) being the cross entropy. Eqn. (8) can be solved for the organ pose Θ( j )* using any gradient-based optimization method. Here we used the steepest descent method [15]. Different

6 y

Centroid (mm)

x 0.06

0.06

0.04

0.04

0.04

0.02

0.02

0.02

Orientation (degree)

0 80

100

120

140

0 120

160

150

180

210

240

0 110

0.06

0.06

0.06

0.04

0.04

0.04

0.02

0.02

0.02

0 -45

-30

-15

0

15

25

0 -30

-15

0

15

30

0 -45

0.06

0.06

0.04

0.04

0.04

0.02

0.02

0.02

0.06

Scaling (%)

z

0.06

0 70

80

90

100

110

0

80

100

120

140

0 40

140

170

200

230

-30

-15

0

15

60

80

100

120

25

Fig. 2. The pose distribution functions of liver. From top to bottom are the centroid, orientation, and scaling respectively; the x, y, and z components are shown from left to right.

organs are localized independently. A multi-resolution strategy is adopted to reduce the computation. In addition, the organ center, orientation, and scaling factor are optimized sequentially for faster convergence and to avoid local optima.

3

Results

The method was applied to abdominal non-contrast CT scans of 12 (6 males and 6 females) patients. Data were collected with a LightSpeed Ultra scanner (GE Healthcare). The slice thickness was 1 mm and the in-slice resolution varied from 0.54 mm to 0.77 mm. In all the 12 images, the five organs (liver, left and right kidneys, spleen, and pancreas) were manually segmented by a medical student and supervised by a radiologist. For validation, 12 experiments were performed with a leave-one-out strategy. In each experiment, one dataset was picked as the subject image, and the remaining 11 datasets served as the training data and were used to build the OPDM and PA. After the organ pose was estimated, the probability atlas was transformed to the subject image using the estimated pose. For quantitative evaluation, the pose errors were calculated by comparing the localized poses with true poses, which were computed from the linear registration as described in Sec. 2.3. The Dice’s coefficient [16] was used to measure the symmetric volume overlap between the estimated organs and manual segmentations. Fig. 2 shows the pose distribution functions, exemplified for the liver, computed using all 12 data sets. Note the large variance of the distributions especially for the yand z-positions and z-scaling. Fig. exemplifies the organ localization results on one dataset using selected axial slices. Error! Reference source not found. presents the

7

Fig. 3. The localization results of five organs on one data set. Images are axial slices from the image volume with the organs identified using different colors. Table 1. The pose estimation errors between the automated localization results and the true poses. In brackets we present the kidney pose errors without considering the outlier cases.

Pose Organ Liver Spleen Left Kidney

Right Kidney Pancreas

mean std mean std mean std mean std mean std

Centroid (mm) x y z 2.42 2.23 4.61 2.38 1.66 2.86 5.28 4.62 3.70 5.63 4.11 3.35 1.06 1.78 9.63 (0.73) (1.34) (2.17) 1.00 1.50 17.7 (0.56) (0.98) (1.73) 1.31 2.69 8.69 (1.08) (1.99) (6.40) 3.24 10.47 1.31 (1.07) (2.23) (7.16) 5.64 6.39 4.17 4.78 7.38 3.68

Orientation (degree) x y z 4.09 2.16 3.50 3.09 1.89 2.61 6.89 5.94 4.93 6.05 4.56 3.72 3.77 4.01 3.90 (3.98) (2.91) (3.48) 4.35 3.72 4.45 (4.74) (2.92) (4.20) 3.13 4.42 3.45 (2.40) (4.20) 3.63) 3.26 3.88 3.92 (2.17) (3.99) (4.06) 6.23 5.05 6.09 6.74 4.62 4.41

x 4.19 3.70 8.31 6.60 6.81 (6.27) 4.41 (4.66) 7.35 (7.00) 3.61 (3.58) 9.20 7.45

Scale (%) y z 5.94 18.8 4.98 6.83 18.8 24.2 8.85 17.2 7.18 4.29 (7.62) (4.72) 5.80 3.27 (6.13) (3.31) 10.04 4.58 (9.93) (4.41) 7.69 3.84 (8.06) (3.99) 10.89 7.79 12.71 7.64

organ pose errors between the organ localization results and the true poses computed according to Section 2.3. Over all, the method performed well for the liver, spleen, left and right kidneys, but less adequate for the pancreas, an organ with higher variability across subjects. In particular for the liver and spleen, the method worked consistently well on all the 12 experiments with small variations. There is a visible correlation between the magnitude of the liver pose errors and the variance of the distributions in Fig. 2. For the left kidney, there were two cases where the kidneys partly overlapped with the nearby spleen (outliers). This was caused by the very similar HU value profiles of spleen and kidneys in non-contrast CT, and their spatial proximity. Left kidneys were well localized in the other 10 cases. For the right kidney, there was one outlier where the method partly converged to the liver for similar reasons. These outlier cases suffered from high location errors in the zdirection of the centroid estimation. Table 1 also shows that the scale estimation has a larger error than the centroid and orientation especially for the liver and the spleen. This is mainly because the scales of these organs present a large variability in the training dataset, as exemplified for the

8 Table 2. The symmetric volume overlap (Dice’s coefficient) after organ localization

mean std

Liver 0.77 0.03

Spleen 0.69 0.06

Left Kidney 0.73 0.21

Right Kidney 0.74 0.14

Pancreas 0.43 0.15

liver scale in z-axis for liver (Fig. 2). We expect these errors to be reduced with a larger training data set. shows the Dice’s coefficient results for the five organs from the leave-one-out experiments. Dice’s coefficient is a common measure of segmentation accuracy with a value of 1.0 meaning that the localized organ completely overlaps with the manual segmentation. Please note that our method localizes abdominal organs and does not address segmentation, hence the inherent inter-patient organ shape variability affects the percentage of overlap. For the left kidney, the two outlier resulted in very low Dice values (0.21 and 0.33), while the other 10 cases had a Dice’s coefficient of 0.82±0.03. The pancreas had a low Dice value in general because of its thin and variable shape. Moreover, the manual segmentation of pancreases was less reliable as it is often difficult to differentiate the pancreas from surrounding tissues in CT data without contrast. Our method was implemented in C++ and run on a computer with 2GB RAM and 3GHz CPU. In our implementation, the localization of liver took about 1 minute, and the localization of each of the other four organs took about 20 seconds. The liver localization takes longer because it has a larger size and thus requires more computation (see Eqn. (6)).

4

Discussion

The proposed organ localization method employed a pose distribution model to limit the searching space, and an MAP framework to estimates the locations, orientations and scales of abdominal organs. The organ localization can be used as an initialization of other medical image analysis tasks in abdominal CAD, e.g., segmentation and registration. The results showed a reasonably good overlap with manual segmentations on liver, spleen, and kidneys, and suggested it can work as an accurate initialization of segmentation or registration tasks. Comparing with other organ localization methods [1,12,13], the proposed method estimates the organ orientations and scales in addition to the locations and thus produces more accurate results. For example, although Yao and Summers’ method [12] estimated organ locations with similar accuracy as the proposed method, their method resulted a much poorer organ localization results because of the lack of knowledge in organ orientations and scales (by visually checking the volume overlaps from Fig. 4 in [12]). The method performed consistently well on the liver and the spleen because they have larger sizes, which helps the MAP optimization to converge correctly. The failed

9

cases of kidneys were mainly due to their smaller sizes and proximity to organs with similar intensity profiles. In addition, the localization of pancreas was more challenging because of its irregular shape and large pose variance across subjects. In this work, different organs were localized independently. In future work, we will further improve the method through modeling inter-organ relations and interactions so that smaller organ localizations can benefit from the localization results of their neighbors. In the future we will use more data to build the models and expect to achieve a more precise pose distribution function and more accurate results. In addition, in noncontrast CT images, the HU values of abdominal soft tissues are similar, which makes it difficult to correct the localization and differentiate organs based on HU values. Better localization results are expected when applying the method to contrastenhanced CT images, and it is our future work. In addition, we will incorporate the proposed method with atlas-based segmentation methods for fully automatic organ segmentation. In conclusion, we developed a novel abdominal multi-organ localization method from non-contrast CT images using an MAP framework. The method computes the organ pose distributions, organ probabilistic atlases, and organ intensity profiles from training data after the abdominal cavity is normalized. It then formulates an MAP problem and solves it for organ poses in given CT images. The method was applied to five organs (liver, left and right kidneys, spleen, and pancreas) with promising results.

Acknowledgements This work was supported by the Intramural Research Program of the National Institutes of Health, Clinical Center. The authors thank Jesse K. Sandberg and Javed Aman for helping with the data analysis.

References 1. Okada, T. et al.: Automatic segmentation of the liver from 3D CT images using probabilistic atlas and multilevel statistical shape model, Academic Radiology 15, pp. 1390-1403, 2008. 2. Song, Y. et al.: Liver segmentation using automatically defined patient specific B-spline surface models, MICCAI 2009 Part II, LNCS 5762, pp. 43-50, 2009. 3. Linguraru, M.G. et al.: Atlas-based automated segmentation of spleen and liver using adaptive enhancement estimation. Med. Phys., Vol. 37(2), pp. 771-783, 2010. 4. Heimann, T., Meinzer, H.-P.: Active shape models for a fully automated 3D segmentation of the liver – an evaluation on clinical data, MICCAI 2006, LNCS 4191, pp. 41-48, 2006. 5. Shimizu, A. et al.: Segmentation of multiple organs in non-contrast 3D abdominal CT images, Int. J. Comp. Assist. Radiol. Surg., 2, pp. 135-142, 2007.

10 6. Linguraru M.G., Summers, R.M.: Multi-organ segmentation in 4D contrast-enhanced abdominal CT, IEEE Int. Symp. Biomed. Imaging (ISBI), pp. 45-48, 2008. 7. Okada, T., et al., Construction of hierarchical multi-organ statistical atlases and their application to multi-organ segmentation from CT images. MICCAI 2008, pp.502-9, 2008. 8. Park, H. Bland, P.H. Meyer, C.R.: Construction of an abdominal probabilistic atlas and its application in segmentation. IEEE Trans Med Imaging, Vol. 22(4), pp. 483-92, 2003. 9. Reyes, M. et al.: Anatomical Variability of Organs via Principal Factor Analysis from the Construction of an Abdominal Probabilistic Atlas. In: IEEE International Symposium on Biomedical Imaging (ISBI), pp. 682-685, 2009. 10. Hu, X. et al.: Independent analysis of four-phase abdominal CT images, Proceedings of MICCAI 2004, LNCS 3217, pp. 916-924, 2004. 11. Seifert, S. et al.: Hierarchical parsing and semantic navigation of full body CT data, Proc. SPIE, vol. 7259, pp. 725902-725902-8, 2009. 12. Yao J., Summers R.M.: Statistical location model for abdominal organ localization. MICCAI 2009, Part II, LNCS 5762, pp. 9-17, 2009. 13. Yao C., et al.: Simultaneous location detection of multi-organ by atlas-guided eigen-organ method in volumetric medical images, Int. J. Comp. Assist. Radiol. Surg. 1, pp.42-45, 2006. 14. Parzen E.: On the estimation of a probability density function and the more, Annals of Math. Stat., 33:1065-1072, 1962. 15. Burden R.L., Faires J.D., Numerical analysis, Brooks Cole, 2004. 16. Dice L.R.: Measures of the amount of ecologic association between species, Ecology 26(3), pp. 297-302, 1945.