Int J CARS (2007) 2:135–142 DOI 10.1007/s11548-007-0135-z
ORIGINAL ARTICLE
Segmentation of multiple organs in non-contrast 3D abdominal CT images Akinobu Shimizu · Rena Ohno · Takaya Ikegami · Hidefumi Kobatake · Shigeru Nawano · Daniel Smutek
Received: 11 April 2007 / Accepted: 31 October 2007 / Published online: 27 November 2007 © CARS 2007
Abstract Objective We propose a simultaneous extraction method for 12 organs from non-contrast three-dimensional abdominal CT images. Materials and methods The proposed method uses an abdominal cavity standardization process and atlas guided segmentation incorporating parameter estimation with the EM algorithm to deal with the large fluctuations in the feature distribution parameters between subjects. Segmentation is then performed using multiple level sets, which minimize the energy function that considers the hierarchy and exclusiveness between organs as well as uniformity of grey values in organs. To assess the performance of the proposed method, ten non-contrast 3D CT volumes were used. Results The accuracy of the feature distribution parameter estimation was slightly improved using the proposed EM method, resulting in better performance of the segmentation process. Nine organs out of twelve were statistically improved compared with the results without the proposed parameter estimation process. The proposed multiple level sets also boosted the performance of the segmentation by 7.2 points on average compared with the atlas guided This study was originally presented at CARS 2006 and supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan. A. Shimizu (B) · R. Ohno · T. Ikegami · H. Kobatake Tokyo University of Agriculture and Technology, Koganei, Tokyo, Japan e-mail:
[email protected] S. Nawano International University of Health and Welfare, Minato-ku, Tokyo, Japan D. Smutek 1st Medical Faculty, Charles University, Prague, Czech Republic
segmentation. Nine out of twelve organs were confirmed to be statistically improved compared with the atlas guided method. Conclusion The proposed method was statistically proved to have better performance in the segmentation of 3D CT volumes. Keywords Computer-assisted diagnosis · Three-dimensional image · Computed tomography · Segmentation of abdominal organs Abbreviations MAP Maximum a posteriori CAD Computer-assisted diagnosis EM Expectation maximization LSM Level set method 3D Three-dimensional IVC Inferior vena cava PV Portal vein Introduction The spatial range and resolution of medical images improved markedly upon the development of multi-slice CT scanners. Now, the entire abdominal area can be scanned with less than a 1-mm-slice interval within a single breath-hold. Although such volume images do aid the diagnosis of multiple organs, a large number of slice images need to be interpreted, leading to lesions occasionally being missed by radiologists. In order to assist radiologists, we have been developing computer-assisted diagnosis (CAD) systems for lesions in the multiple abdominal organs [1]. Accurate segmentation of the multiple abdominal organs is very important for such CAD systems from the viewpoint of diagnostic accuracy and computational cost. To this end,
123
136
Int J CARS (2007) 2:135–142
Fig. 1 Organs to be extracted. Boundaries of different organs are shown in different colours. To evaluate the performance of the proposed method, all the organ boundaries were manually defined by the authors and approved by a medical doctor
extraction algorithms that can process the multiple abdominal organs have been proposed. For example, Park et al. [2] proposed an atlas-based segmentation of abdominal organs using a combination of the probabilistic atlas and the maximum a posteriori (MAP) approach, where they demonstrated successful integration of the atlas information into the segmentation. However, segmentation results have irregular shapes due to the large variation in shape and grey value of organs among different human bodies. Recently, there have also been significant efforts to use the spatial relationship between regions to boost the performance of multiple objects segmentation [3], constraints by neighbouring regions [4], and interactions between neighbouring organs [5], all of which are very effective in extracting the multiple structures in the brain and pelvis in MRI. Although these methods can be applied to abdominal organs, less attention has been paid to actually segmenting the multiple abdominal organs in non-contrast CT images, as there are several important differences between brain and abdominal images. Namely, abdominal organs are mainly surrounded by soft tissues and show large variations in their positions and shapes, which is a different case from the brain within the solid skull. Consequently, it has become necessary to develop a new algorithm that can reduce these variations and consider the spatial relationship between the abdominal organs to obtain fine segmentation results. Park et al. proposed the standardisation process for four organs—liver, kidneys and spinal cord. They used the thin-plate spline based method to align the input data with the target reference geometry. Such nonrigid registration is effective in reducing variation not only in position but also in size and shape of organs. However, the non-linear registration technique may cause unfavourable distortion in the region where the non-continuous transformation occurs, such as sliding between neighbouring organs, which makes the segmentation task more difficult. In this study, we focus on 12 organs with such a non-continuous deformation. To prevent excessive distortion, we employed linear transformation to standardize the abdominal cavity. The proposed standardization method is based on automati-
123
cally detected landmarks. Subsequently, it performs an atlas guided segmentation algorithm incorporating an expectation maximisation (EM) algorithm for the extended-mixture Gaussian distribution in which the mixture ratio is defined at each voxel. We also propose a multiple level set method (LSM), which minimises the energy function that considers the hierarchy and exclusiveness between organs as well as the uniformity of grey values in organs. The LSM is implemented in a clustered computer environment using message passing interface to reduce the computational cost. We report here the experimental results of applying the algorithm to 12 organs using actual non-contrast 3D CT images. Materials and methods Materials Ten abdominal non-contrast 3D CT images taken from ten males were used to assess the performance of the proposed method, where the average number of slices per case was 184, the slice interval was 1 mm and the spatial resolution in the axial slices was 0.625 mm. The images were converted to isotropic voxels and image size was reduced to one-third for the following experiments. Target organs were the liver, spleen, left and right kidneys, heart, gallbladder, pancreas, portal vein (PV) between the spleen and liver, stomach wall, oesophagus, abdominal aorta and inferior vena cava (IVC) (Fig. 1). Segmentation of multiple organs from abdominal CT images Figure 2 shows a flowchart of the proposed segmentation process consisting of standardisation of the abdominal cavity, rough extraction and fine segmentation. Standardization of abdominal cavity Although many papers have described standardization of the brain [6], not much attention has been paid to the
Int J CARS (2007) 2:135–142
137
Probabilistic atlas
Feature database
3D CT image
Standardization of abdominal cavity Atlas guided rough extraction Fine segmentation based on multiple level set method Segmented regions
Fig. 2 Outline of the proposed algorithm, which has three steps: standardization of the abdominal cavity, rough segmentation based on a digital atlas of human anatomy, and fine segmentation using the multiple level sets
the changes in the lung area in axial slices, z coordinates of the right lung base, and directional difference in CT values along the body axis around the right lung base. The median of the estimated positions by these three detectors is defined as the top of the diaphragm. Next, the remaining four landmarks are detected independently using the template matching technique, which maximises the correlation coefficient between the input image and the average templates of the kidneys, liver, and spleen. It exhaustively searches within the restricted area, determined statistically from the training dataset, and then the positions with maximum correlation are defined. Finally, we apply an affine transformation, which minimizes the difference between the detected landmarks of the input image and those in the reference image. The minimization process is performed using Powell’s method [7].
Atlas guided segmentation standardization of the abdominal cavity. Since abdominal organs are mainly surrounded by soft tissues, there are large variations in the positions and shapes of organs. Park et al. proposed a standardisation process based on the non-linear registration technique, which is effective in reducing variation in the location of organs as well as in their shape and size. However, such a flexible method might cause excessive distortion in the region where non-continuous transformation occurs and result in decreasing the performance of the segmentation. Therefore, we used an affine transformation to normalize the abdominal cavity. The abdominal cavity includes several landmark points, which are relatively easy to detect and help to normalize the abdominal cavity, such as the tops of the kidneys. Therefore, we have focused on a landmark-based standardization process and employed five landmarks that are relatively easy to detect and are widely distributed in the abdominal cavity. The landmarks are the top voxels of the diaphragm and kidneys, as well as the lower end voxels of the liver and spleen (Fig. 3). First, the top of the diaphragm is determined by the detector ensemble, each part of which estimates the position using
Inputs for this atlas-based method are 3D abdominal CT images and a digital atlas of human anatomy, which consists of a probabilistic atlas [2] and feature database including feature distribution parameters, such as average and covariance. Here the probabilistic atlas is computed by applying the same standardisation to each manually segmented image. Then, for each voxel in the reference volume, the fractional percentage of data that has a label at that reference voxel location is calculated. Figure 4 shows volume rendering of the probabilistic atlas of liver generated from manually segmented liver regions in 3D CT images obtained from ten patients. In the following equations, the label space is denoted L and the observed data (= feature vector) is denoted as V . Sample realizations of L and V are represented as lower case letters. The atlas guided segmentation is realized using the maximum a posteriori (MAP) algorithm. The purpose of this algorithm is to estimate the label L that best explains the observation V by maximizing the a posteriori probability given V . lˆ = arg max P(L = l|V = v)
(1)
l
Fig. 3 Five landmarks used for standardisation of the abdominal cavity
123
138
Int J CARS (2007) 2:135–142
Fig. 4 Examples of input images. Left CT images, right probabilistic atlas of liver
Using Bayes’ formula, the a posteriori probability p(l|v) is replaced as follows: P(L = l|V = v) =
p(v|l) p(l) p(v)
(2)
The elements of V were (x, y, z) coordinates at the voxel of interest and the CT values of the voxel; p(l) shows a priori probability of organ l, which is derived from the probabilistic atlas. In order to cope with the changes of the distribution parameters of the feature vector V between subjects, our method estimates the parameters by the EM algorithm for the extended mixture Gaussian distributions as follows: 1 αl (n)N(v; µl , l ) p(v; θ ) = N M
N
(3)
l=1 n=1
where M and N are the number of organs and voxels in the image, respectively. N (v; µl , l ) is a Gaussian distribution with average µl and covariance l for organ l. Compared with the conventional mixture model, the mixture ratio αl (n) is extended to be defined at each voxel, which is consistent with the fact that each organ occupies a limited area. Moreover, such extension allows us to use the prior probability of organs in the digital atlas of human anatomy, resulting in higher accuracy in parameter estimation than in conventional mixture distribution, which will be demonstrated in the next section. We have derived the update equation of the mixture ratio αl (n) by introducing the Lagrange multiplier λ with the M N constraint l=1 n=1 αl (n) = N resulting in (t+1) (4) (n) = p l|xn , θ (t) αl Note that there are no changes from the update equations for µl and l of the conventional mixture distribution [8].
123
Finally, the proposed method applies morphological operations or opening and closing with the structure element whose size was optimized experimentally to eliminate unnecessary protrusions, holes and cavities. We performed the experiment to determine the optimal size for each organ. We applied the structure elements with different sizes to each organ extracted by the atlas-guided segmentation and selected the size when the accuracy (Eq. (9) in the next section) of the extracted regions from the training dataset was maximal. Note that the optimal size depends on organ but not patient. So it is not necessary to change the size when you apply the morphological operations to different patients.
Fine segmentation based on the multiple level set method Although the above process can roughly delineate the organs’ borders, details of the boundaries are missed. Thus, the proposed method sophisticates the boundaries using the level set method (LSM). The introduction of constraints from neighbouring organs is known to be effective in obtaining fine segmentation [3–5]. Moreover, variation in the CT value of the exactly extracted region is usually lower than that of incorrect regions, which include several organs of different CT values. Consequently, we have developed a new method for extracting abdominal organs, which minimises the energy function that considers the hierarchy and exclusiveness between organs and uniformity of grey values in the organs. Energy relating to the hierarchy of the organs is defined as overlapping volume with outside regions of the abdominal cavity (Fig. 5) extracted with an active cylinder model [9]. Energy associated with exclusiveness is designed based on overlapping volume with neighbouring organs. Deviation in grey values I (x, y, z) is also introduced into the energy function by referring [10]. The proposed energy function
Int J CARS (2007) 2:135–142
139
c=
1 1 + |∇(Gσ ∗ I)|
where G σ is a Gaussian operator with standard deviation σ .
Results
Fig. 5 Estimated outside region of abdominal cavity using active cylinder model [9]
expressed using regularized delta and Heaviside functions follows: ⎡ ⎡ M M ⎣ ωi j ⎣αi j (1 − Hε (φi )) E= j=1, j=i
i=1
⎤
error = |estimated value–true value|
×(1 − Hε (φ j ))d xd ydz ⎦ +βi
The number of Gaussian distributions assumed in the atlasguided segmentation is 13, including 12 organs and another tissue class corresponding to the non-target organ, such as bowel. Note that the non-organ regions, such as fat and air, as well as the lungs and bones are extracted in advance using conventional image processing, or binarisation followed by morphological processing. The LSM was implemented in a clustered computer environment that included 16 CPUs, each of which had a 2.8 GHz clock speed with 4 GB memory and was assigned to one of the organs to be extracted. The validation test was done by the re-substitution method where the dataset for learning is the same as that for the test. Estimation accuracy of the feature distribution parameters (average µl and covariance l ) in the atlas guided segmentation process was evaluated by
⎤ (I(x, y, z) − µi )2 (1 − Hε (φi ))d xd ydz ⎦ (5) σi2
where ωi j , αi j and βi are weights and φi is a level set function of organ i. shows an image domain, and µi and i are the estimated average and variance of CT values. We deduced the associated Euler–Lagrange equation, which minimizes the energy function with respect to φi , while keeping µi and i fixed. Parameterizing the descent direction by artificial time t (>0), the evolution equations in φi (i = 1, 2, . . . , N ) are ⎡ M
∂φi = δε (φi ) ⎣ ωi j αi j 1 − Hε φ j ∂t j=1, j=i (I(x, y, z) − µi )2 +βi (6) σi2 In this study, we linearly combine the function with the following conventional evolution equation [11] in order to consider the edges and the curvature in the process of evolution: ∂φi (7) = δε (φi )[c(κ + V0 )|∇φ| + ∇c · ∇φ] ∂t where κ is a mean curvature of a surface of organ and V0 is a constant value that makes a level set function go up or down with a constant speed.
(8)
where the true values were defined using the true regions of organs manually segmented by authors and approved by a medical doctor. Note that the µl and l in Eq. (3) depend on organs but not on the location n. So the error of Eq. (8) was evaluated at each organ. Figure 6 demonstrates that estimation accuracy was improved by the proposed estimation process denoted as EM-ex: EM algorithm for extended mixture Gaussian distribution. In the case of liver, both the conventional and proposed estimation processes gave good results. The errors were less than two. For other organs, however, the errors of the proposed process were smaller than those of the conventional one. The errors were 32 points lower on average than the errors with the EM algorithm for conventional mixture distribution, which significantly failed to estimate the parameters of the small or elongated organs, such as the pancreas PV and IVC. The number of points in the above sentence was computed by dividing the difference between the errors of the two processes by the error of the conventional one. Compared to the errors of the parameters stored in the feature database (denoted as DB in the figure), the number of errors in the proposed method was 23 points fewer on average. Figure 7 presents examples of the extracted boundaries by MAP segmentation, which uses the estimated parameters. This figure shows that, although, most of the organs were roughly extracted, underestimation and overestimation still occurred compared with true regions when focusing on regions near boundaries.
123
Int J CARS (2007) 2:135–142 20 18 16 14 12 10 8 6 4 2 0
DB EM
C IV
ta A or
s ag u ph
O es o
St om ac h
PV
er Pa nc re as
al lb la dd
G
H ea rt
ey
ey
K id n R.
id n L. K
Sp
le en
EM-ex
Li ve r
Error
140
Fig. 6 Mean and standard deviation of errors between the estimated average and true average. DB indicates errors between parameters in the database and true parameters. EM and EM-ex show results of the EM algorithm for the conventional mixture model and the extended mixture model, respectively
Fig. 7 Examples of extracted boundaries by atlas guided segmentation
Fig. 8 Left and middle figures show results without/with the new energy term, respectively. Please note that only the boundaries of liver and right kidney are indicated for sake of simplicity. The right figure depicts a 3D view of the extracted regions of all organs by the proposed method
We evaluated the performance of MAP segmentation based on the degree of correspondence between the results and the corresponding manually segmented true regions, as defined by degree of correspondence (%) =
#(extracted region ∩ true region) × 100 #(extracted region ∪ true region)
123
(9)
where #(region) is a function that computes the number of voxels in the region. The degree of correspondence ranged from 0 (completely different) to 100 (completely overlapped). Using this validation test, we found that the degree of correspondence of the proposed method was four points better on average than the segmentation results using the parameters in the database. It was still higher by 2 points on average than the segmentation results with the EM algorithm for conventional mixture distribution.
Int J CARS (2007) 2:135–142
141
Degree of correspondence
1 0.9
MAP
0.8
conventional energy
0.7
proposed energy
0.6 0.5 0.4 0.3 0.2 0.1
C IV
ta or
op
A
ha gu s
h O
es
PV
s ea
m ac St o
la lb G
al
cr
er dd
ea H
Pa n
rt
y K
id
ne
y R.
id
ne
n ee L. K
Sp l
Li v
er
0
Fig. 9 Performance of the segmentation process for each organ. MAP gives initial boundary of the level set method. This figure shows that the proposed energy function outperforms the conventional energy function
Fig. 10 Example of the segmentation by atlas guided segmentation. Left figure shows the true boundaries manually defined by authors. Middle one corresponds to the result based on the conventional parameter estimation algorithm. The right one shows the result using the proposed algorithm
The subsequent LSM refined the boundaries of the organs. Figure 8 shows resultant examples without/with the proposed evolution equation. The red lines show the resultant boundaries of the liver, while the orange lines correspond to those of the right kidney. The boundaries in the middle figure successfully demonstrate the improvement by the proposed LSM. The degree of correspondence computed from Eq. (9) was improved by six points for liver and two points for the right kidney. Figure 9 presents the degree of correspondence for each organ. Based on this figure, the degree of correspondence was increased by 3.7 points on average compared with MAP based segmentation and 3.5 points better by introducing the new energy term. However, the boundaries of elongated organs with blurred boundaries, such as the pancreas, oesophagus, IVC and PV were poorly defined. The average degree of correspondence of these four organs with the stomach wall was 32.5%. In contrast, the average performance of the remaining seven organs (liver, spleen, left and right kidneys, heart, gallbladder and abdominal aorta) was 78.8%, which is acceptable for a subsequent diagnosis process.
Discussion In this paper, we have presented a novel algorithm for the simultaneous extraction of 12 organs from non-contrast 3D abdominal CT images. The main contributions of this paper are the extension of the MAP based segmentation algorithm to be applicable to 12 organs in non-contrast CT images and the proposal of EM algorithm for extended mixture Gaussian distribution. Moreover the concatenation of MAP based segmentation to LSM with mutual interactions between neighbouring level sets is another important proposal. With the proposed EM algorithm, unfavourable fluctuation of the features between subjects was reduced. A typical example of the atlas-guided segmentation is presented in Fig. 10 to show the advantage of using the proposed parameter estimation algorithm. The regions denoted by yellow arrows suggested to us that the segmentation algorithm based on the proposed EM algorithm outperformed the conventional EM algorithm. The accuracy improvement in the parameter estimation process resulted in a high accuracy estimation of the regions. To statistically assess the difference in performance (=“degree of correspondence” defined
123
142
Int J CARS (2007) 2:135–142
Fig. 11 Example of the segmentation by LSM. Left figure shows the true boundaries. Middle one corresponds to the result based on the conventional energy, and the right one shows the result using the proposed energy
by Eq. (9)) statistically, we carried out a single-tailed paired t-test. The null hypothesis was no improvement in performance. The significance level was 0.05. When comparing the segmentation results using the parameters in the database with that of the proposed algorithm, we found that the improvement was statistically significant for nine organs and no difference was observed for three organs. On the contrary, when using conventional algorithm, there were two organs whose performance was significantly decreased, while eight organs were improved. The same statistical test was carried out to evaluate the improvement by the proposed LSM. Comparing the degree of correspondence in Fig. 9 of the LSM using the conventional energy with those of the atlas-guided segmentation, accuracy in the segmentation of five organs was improved, but performance of oesophagus segmentation was decreased significantly. In contrast, the LSM based on the new energy increased the performance in segmentation of nine organs, and no deterioration was observed. The statistical difference between the results using the conventional energy and those of the proposed energy was found in seven organs. An example is presented in Fig. 11 where mutual interaction improved the accuracy in segmentation of organs as indicated by the yellow arrows. However, such interaction occasionally caused error in segmentation. The failure to extract the abdominal cavity’s boundary led to the error in segmentation of the heart as indicated by the red arrow. As shown in Fig. 9, the proposed method was not sufficiently accurate in the segmentation of the small and elongated organs, such as the oesophagus, PV and pancreas, even though it is difficult for radiologists to define these boundaries on non-contrast CT images. When segmenting such structures, a priori knowledge of the shapes and spatial relationships between organs is very helpful [3–5,12]. In the future, a novel a posteriori probability that takes into account the probabilities of neighbouring organs will be introduced
123
into the atlas-guided segmentation process. We also plan to develop a new energy term to incorporate the shape a priori information to recover the 3D shapes. Testing the performance, using a large-scale database and extension of the method to be applicable to the whole body are other exciting and challenging future projects. References 1. Kobatake H (2007) Future CAD in multi-dimensional medical images - project on multi-organ, multi-disease CAD system - Comput Med Imaging Graph 31(4–5):258–266 2. Park H, Bland PH, Meyer CR (2003) Construction of an abdominal probabilistic atlas its application in segmentation. IEEE Trans Med Imaging 22(4):483–492 3. Yang J, Staib L, Duncan J (2004) Neighbor-constrained segmentation with level set based 3-D deformable models. IEEE Trans Med Imaging 23(8):940–948 4. Tsai A, Wells W, Tempany C et al (2003) Coupled multi-shape model and mutual information for medical image segmentation. Inf Process Med Imaging 18:185–197 5. Yan P, Shen W, Kassim AA et al (2005) Segmentation of neighboring organs in medical image with model competition. In: Proceedings of MICCAI 2005, LNCS3749, pp 270–277 6. Maintz JBA, Viergever MA (1998) A survey of medical image registration. Med Image Anal 2:1–36 7. Press WH, Flannery BP, Teukolsky SA et al (1995) Numerical recipes in C, 2nd edn. Cambridge University Press, Kluwer, London, Dordrecht, pp 412–420 8. Ghahramami Z, Jordan M (1995) Learning from incomplete data. Technical Report AI Lab Memo no. 1509, CBCL Paper no. 108, MIT AI Lab 9. Okumura T, Yamamoto S, Matsumoto M et al (1998) The lung region extraction in the chest CT images by the active cylinder model. Med Imaging Technol 16(1):61–71 10. Chan TF, Vese LF (2001) Active contours without edges. IEEE Trans Image Process 10(2):266–277 11. Caselles V, Kimmel R, Sapiro G (1997) Geodesic active contours. Int J Comput Vis 22(1):61–79 12. Yao C, Wada T, Shimizu A et al (2006) Probabilistic atlas-guided Eigen-organ method for simultaneous bounding box estimation of multiple organs in volumetric CT images. Med Imaging Technol 24(3):191–200
