This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging IEEE TRANSACTIONS ON MEDICAL IMAGING
1
Automated Segmentation of Breast in 3D MR Images Using a Robust Atlas Farzad Khalvati, Cristina Gallego Ortiz, Sharmila Balasingham, Anne L. Martel
Abstract—This paper presents a robust atlas-based segmentation (ABS) algorithm for segmentation of the breast boundary in 3D MR images. The proposed algorithm combines the wellknown methodologies of ABS namely probabilistic atlas and atlas selection approaches into a single framework where two configurations are realized. The algorithm uses phase congruency maps to create an atlas which is robust to intensity variations. This allows an atlas derived from images acquired with one MR imaging sequence to be used to segment images acquired with a different MR imaging sequence and eliminates the need for intensity-based registration. Images acquired using a Dixon sequence were used to create an atlas which was used to segment both Dixon images (intra-sequence) and T1-weighted images (inter-sequence). In both cases, highly accurate results were achieved with the median Dice Similarity Coefficient (DSC) values of 94% ± 4% and 87 ± 6.5%, respectively. Index Terms—Breast MRI, Image Segmentation, Atlas-Based Segmentation, Image Registration.
I. I NTRODUCTION
M
AGNETIC resonance imaging (MRI) of breast is an increasingly common approach for monitoring and detection of breast cancer mainly due to higher sensitivity and no ionizing radiation compared to conventional x-ray mammography. MR imaging is now also recommended for screening women who are known to be at a higher risk of breast cancer [1]. The segmentation of breast tissue in MR images has a number of important applications. The removal of unwanted structures such as the chest wall and the heart makes it possible to improve the quality of 3D visualisation. In order to use MRI to assess breast density, a significant risk factor and an important biomarker in determining the possibility of breast cancer [2], it is essential to determine the total volume of the breast. In Computer-Aided Diagnosis (CAD), breast segmentation allows the software to ignore tissue outside of the breast, reducing the false positive rate. Finally, in order to plan therapeutic interventions such as radiotherapy using MRI, accurate segmentations of the target organs will be required [3]. Atlas-based segmentation (ABS) is a well established and widely used technique for extracting contours from medical images of different organs (e.g., prostate) [4], [5], [6], [7], [8], [9]. In this method, processed images are stored in a Copyright (c) 2014 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to [email protected]. F. Khalvati is with the Sunnybrook Research Institute, Toronto, Ontario, Canada. e-mail: [email protected] C. Gallego Ortiz, S. Balasingham, and A. L. Martel are with the Sunnybrook Research Institute, Toronto, Ontario, Canada.
database called a set of atlases. Each atlas consists of an image and its gold standard (i.e., manual segmentation or label). The atlas images are first registered to the target image and then the resulting transformations are applied to the labels. Finally, the deformed labels are combined to achieve the final segmentation result. In general, there are three approaches to design an ABS algorithm: multi-atlas, atlas selection, and probabilistic atlas. The multi-atlas approach registers a limited number of atlas images to the target image to create multiple (deformed) labels, which are then combined using a fusion algorithm to generate the final result. Different fusion algorithms for multi-atlas approaches have been proposed in the literature. Majority voting is the simplest approach where each pixel is assigned to the label which occurs with the highest frequency. More complex methods attempt to weigh the segmentation results from each atlas based on their estimated accuracy [5]. Langerak et al. [9] proposed a multi-atlas segmentation method for prostate MR images based on an iterative label fusion approach. First, all atlas images are registered to the target image to obtain a set of deformed labels. The deformed labels are fused together using a weighted majority voting method and then, the overlap of each label is calculated against the fused label. Labels with low overlaps are discarded and a new fused label is recalculated. As another label fusion method, the STAPLE algorithm [10] uses an expectation-maximization approach in which the unknown true segmentation and the performance of individual atlas segmentations are estimated simultaneously. The atlas selection approach proposed in [5] uses both STAPLE and majority voting as the post-processing stage to fuse the generated labels. Another approach is to assume that the performance of each segmentation is related to the similarity between the registered atlas image and the target image. In [11], it was found that weighting labels using a local image similarity measure given by the mean squared difference between target and registered atlas image gave better results than majority voting or STAPLE. Several atlas selection methods have been proposed which try to find the subset of atlas images that best represent a specific target image in order to improve accuracy. The mutual information (MI) between the atlas and the target image after registration has been used as a figure of merit; for example in [12] the target images are subdivided into regions and for each region the atlas image with the best local MI is selected, and in [5] the subset of atlas images that has an MI greater than a pre-set threshold level is used. Atlas selection approaches that require all atlas images to be
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 2
IEEE TRANSACTIONS ON MEDICAL IMAGING
registered to the target image are computationally prohibitive and limit the number of atlases that can be used. To reduce the computational cost of atlas selection, Aljabar et al. [13] proposed to register all atlas images and the target image to an arbitrary reference image and thus, the best-match atlas is selected by comparing the registered target image with the registered atlas images. In other words, the atlas registration can be performed off-line during the atlas construction phase reducing the computational time. In a recent work, Langerak et al. [14] proposed a method that clusters the atlas images using the pairwise registration results to group similar atlases together. For each cluster, an exemplar atlas is selected and these are then registered to the target image. Once the best exemplar is found, all of the atlas images belonging to that cluster are registered to the target image and the resulting labels are combined using a simple majority voting. Although the process of clustering the atlases off-line is extremely computationally expensive, the number of registrations required for each target image is somewhat reduced. The probabilistic atlas approaches are usually based on either sequential models or generative models. In sequential models [4], [7], the probabilistic atlas is generated by aligning (and averaging) multiple images and their labels in a dataset using a registration algorithm. In generative models [15], [16], [17], [18], the tissue classification and image deformation are performed simultaneously using an objective function optimized by deforming the tissue probability maps. The main difference between the two classes is that in the former, the images are first registered and then combined to generate a probabilistic atlas. Next, the resulting probabilistic atlas is registered to the target image and consequently, the probabilistic label is deformed via the registration transformation to generate the result for the target image. In contrast, in the latter model, the target image is segmented (classified) directly as a result of deforming the tissue maps of the image via optimizing the objective function. There are different variations of sequential models of probabilistic atlas algorithms in the literature. For instance, Martin et al. [4] proposed to create a probabilistic atlas for prostate by registering images to a manually picked reference image. Next, a mean image is created by averaging the registered images. Each image is registered against the mean image to produce a deformed label of the image. The deformed labels are then averaged to generate a probability map of the labels. To segment a target image, the mean image is registered to it and the probability map of the labels is deformed using the registration transformation. Dowling et al. [7] also proposed a similar approach. To build a probabilistic atlas using a sequential model, a reference image(s) which has already been manually segmented and selected plays a crucial role. In practice, however, generating a reference image is not straight forward; even with careful atlas selection, the selected image may not represent the population well enough and hence, the registration may not produce acceptable results. Balci et al. [19] extended the groupwise registration algorithm proposed in [20] to be an alternative to pair-wise registration when designing a sequential probabilistic atlas for segmenting brain MR images.
Groupwise registration aligns a set of images to a common reference image by generating a set of transformations that map the reference image to each of the images in the group. This makes the sequential probabilistic atlas more robust and unbiased to a particular reference image. Ashburner and Friston [16] proposed a generative model of probabilistic atlas which unifies tissue classification, bias correction, and template registration for brain MR images. The objective function is based on a mixture of Gaussian (MOG) models where the tissue classes are represented by Gaussian probabilities and the bias, which is a smooth and spatially varying artifact that corrupts the image, is taken into account by extra parameters in MOG model. In this method, a template (spatial priors) is used in the process of optimizing the objective function where the deformation (registration) is updated during the process. Recently, Iglesias et al. [17] proposed a framework for multi-atlas segmentation of brain in MR images. The proposed algorithm is an extension of the work presented in [16] where in contrast to [16], it builds multiple probabilistic atlases. Although the ABS algorithms have been applied extensively to different organs such as brain and prostate, there have been limited work on breast segmentation using ABS in the literature. Gubern et al. [21] proposed a probabilistic atlas algorithm to segment different structures in breast MR (T1-weighted) images. In this method, the probabilistic atlas is created using pair-wise registration of each atlas with a reference image. The segmentation is done in a Bayesian framework where for a given target image voxel, the segmentation is estimated by searching for the most probable label [21]. The generative models of probabilistic ABS algorithms are usually computationally expensive since they integrate multiple optimizations [16]. The main drawback of the sequential models of probabilistic ABS algorithms is that they usually depend on the image intensity, rather than image key characteristics. As a result, if intensity variations exist between the target image and atlas image, which is common when the two images are from different scanners, the ABS algorithms may perform poorly resulting in inaccurate segmentation results. Moreover, in the conventional ABS algorithms, the atlas and target images must be from the same image sequence (e.g., MRI T1-weighted) in order to obtain reasonable segmentation results. One possible approach to make the sequential models of probabilistic ABS intensity invariant is to use a registration algorithm that is robust to intensity variations. Heinrich et al. [22] proposed an intensity-robust feature for image registration based on modality independent neighbourhood descriptor, which can be used for nonrigid registration across different modalities and different sequences of one modality (e.g., MRI). In this paper, we propose an improved ABS algorithm which is novel in that it combines the atlas selection method with a probabilistic atlas generated using a sequential model. The atlas images are grouped into classes based on their pairwise similarities. Each class of atlases is then used to create a probabilistic atlas using a groupwise registration. The proposed algorithm is presented in two configurations
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
for the intra-sequence and inter-sequence settings where in the former, the atlas and target images are from the same MR imaging sequence and in the latter, they are from different MR imaging sequences. In order to tackle the intensity-dependence deficiency of conventional ABS algorithms, for the intersequence setting, we use phase congruency maps, which are robust to intensity variations, to build the probabilistic atlas1 . These intensity-invariant maps are also used to register the atlas and target images. In order to refine the atlas results, the actual edge information from the target image is also incorporated directly into the segmentation pipeline by using an active contour model. To reduce the computational time, our proposed approach performs the atlas selection in pre-nonrigid registration stage and hence, to segment a target image, it only requires one nonrigid registration operation. II. M ATERIALS The proposed algorithm in this paper has been configured for two settings: intra-sequence and inter-sequence. In the intra-sequence setting, the atlas and test (target) images2 are from the same MR imaging protocol namely Dixon images (section II-A). In the inter-sequence setting, the atlas and test images are from different MR imaging protocols namely Dixon images and MR T1-weighted images, respectively (section II-B). In the following sections, the datasets used in this study are described in more detail. A. Intra-Sequence Materials For the intra-sequence setting, the atlas and test data, used from a previous study [23], contained 400 breast MRI datasets (94×94×44 pixels, wtih a resolution of 2.56×2.56×2.56 mm) acquired using a modified Dixon sequence. The population consisted of two age groups; 320 women aged 15-30 and 80 women aged 40-60 who were healthy and recruited between December 2003 and December 2007. The fat and water images of breast were obtained on a 1.5T Signa Cvi MR system (GE, Waukesha, WI) using a dedicated coil for MR breast (Medical Advances, Milwaukee, WI). The subjects were scanned in the prone position and the images were acquired in the sagittal plane with a slightly modified version of the GE FSE Dixon sequence. The images used in this paper were of zero phase shift. Out of 400 MRI datasets, 350 datasets were selected randomly to create the atlas and 50 datasets were set aside and used to evaluate the segmentation performance. A k-fold cross-validation was performed with k = 5 and the results were averaged. The images were contoured by an expert user by manual initialization of an active contour technique [24] where the segmentation results were corrected by the user to create ground-truth results. B. Inter-Sequence Materials For the inter-sequence setting, the data used to create the atlas was the same as the one used for the intra-sequence 1 In the remaining of the paper, probabilistic atlas refers to a sequential model of the probabilistic atlas approach. 2 In the remaining of the paper, test image and target image are used interchangeably.
3
setting as described in the previous section, which contained 350 breast MRI datasets (94 × 94 × 44 pixels) of Dixon imaging sequence. The test data comprised MR T1-weighted images acquired in 2013 using 3 different MRI scanners at the Sunnybrook Health Sciences Centre, Toronto, Ontario, Canada, for 9 healthy volunteer subjects. The age of the subjects ranged from 33 to 60 years and their BMIs ranged between 20.1 to 30.1 kg/m2 . All data was obtained under the local institutional research ethics board. For 8 subjects, both breasts were scanned while for one subject, only the right breast was scanned. The subjects were scanned in the prone position and the images were acquired in the sagittal plane. In order to evaluate the performance of the algorithm on different scanners, three scanners of both 1.5T and 3T magnets were used with different image dimensions. Table I summarizes the scanners and acquired data used to evaluate the performance of the proposed algorithm. The ground-truth results were created by manual segmentation of the images using TurtleSeg [25] (Interactive 3D Image Segmentation) software. III. P ROPOSED S EGMENTATION A LGORITHM FOR I NTRA -S EQUENCE S ETTING Our ABS algorithm combines a sequential model of probabilistic atlas and atlas selection approaches to offer an efficient ABS method. The goal is to increase the accuracy of the results by diversifying the atlas. In this section, a configuration of the algorithm for the intra-sequence setting is presented where the atlas and target images are from the same MR image sequence. A. Intra-Sequence Setting: Atlas Construction Atlas construction consists of two phases: clustering function and the probabilistic atlas building. 1) Clustering Function: Let F (x) be a target 3D image (or a Fixed image) and A = {a1 , a2 , ..., an } be a set of atlases where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. The first step is to cluster the atlases into N 3 separate classes so that similar atlases are grouped into one class. This requires a similarity measure that calculates the match between any given pair of atlas images in A, Sim(Ii (x), Ij (x)). This builds a (dis)similarity matrix of atlas images where the matrix represents pairwise dissimilarities between all atlas images. For the similarity measure, we use the cross-correlation coefficient4 . Sim(Ii (x), Ij (x)) = Corr(Ii (x), Ij (x))
(1)
For n atlas images, the similarity function Sim creates (dis)similarity values, which are fed to a multidimensional scaling (MDS) algorithm [26] to create a 2D distance map of all images5 . Next, the K-means clustering algorithm [27] is used to cluster images into N classes based n(n−1) 2
3 As it will be discussed in section III-C, for a given training dataset, the optimal value for N can be determined by running experiments across different values of N . 4 Mutual information may also be used as the similarity measure. 5 It is assumed that the atlas images are roughly aligned. Otherwise, an affine registration must be used to align the images before clustering.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 4
IEEE TRANSACTIONS ON MEDICAL IMAGING
TABLE I D ESCRIPTION OF THE BREAST MR T1- WEIGHTED IMAGES USED FOR INTER - SEQUENCE SETTING Total studies 8 4 2 2 1
Scanner GE MR7250 3T, GE MR7250 3T, Phillips Phillips
Dimensions (pixels)
GE Signa CV/i 1.5T, 8 channel coil Sentinelle 16 channel Vanguard Breast MRI tabletop coil Sentinelle 16 channel Vanguard Breast MRI tabletop coil Achieva 3T, Philips 7 channel Sense Breast coil Achieva 3T, Philips 7 channel Sense Breast coil
on the distance map where each class of atlases Acj has mj atlases where j ∈ [1, N ]. Acj = {a1 , a2 , ..., amj }
(2)
The union and intersection of atlas classes will give the set of all atlases A and an empty set, respectively. 2) Probabilistic Atlas: At this point, each class consists of several atlas images with the corresponding labels. The next step in atlas construction is to create a probabilistic atlas for each class of atlases. To create the probabilistic atlas, we use a groupwise registration technique which aligns a set of images to a common reference image. Balci et al. [19] proposed a framework for groupwise registration of images which uses the stack entropy cost function and a multi-resolution B-spline nonrigid deformation to generate the set of image transformations. The motivation behind using stack entropy as a cost function was that by aligning the images accurately, the intensity values of pixels in the corresponding locations in the stacked images should not vary significantly, which means that the stack entropy should be low [19]. The proposed approach uses a combination of global and local transformations where the former is an affine transformation and the latter is a nonrigid deformation based on B-splines [28]. In order to create the probabilistic atlas for each class, first, the groupwise registration is applied to atlas images in each class to generate the probabilistic image for the class (Equation 3). The transformations are applied to labels of each class to generate the probabilistic label of the class (Equation 4). Ijp (x) =
mj 1 X k Ij (x) ◦ TIjk (x)→IjR (x) mj
(3)
mj 1 X k Lj (x) ◦ TIjk (x)→IjR (x) mj
(4)
k
× × × × ×
0.3516 0.3516 0.4688 0.9375 0.5804
× × × × ×
2 2 2 2 2
Lpj (x)
=
1 mj
k=1 m Pj
k=1
Lkj (x) ◦ TIjk (x)→IjR (x)
B. Intra-Sequence Setting: Atlas Matching In the atlas matching phase, the target image F (x) is first aligned (via affine registration [29]) and compared to the probabilistic atlas images to find the best-match atlas image p ImaxSim (x). p f ine ImaxSim (x) = arg maxj Corr(F (x), Ijp (x) ◦ TIaf ) p j (x)→F (x) (5) Once the best-match probabilistic atlas image is found, it is registered (via nonrigid registration [29]) against the target image to obtain the registration transformation. The obtained transformation is applied to the best-match atlas label to generate the segmentation result LA F (x). p nonrigid LA F (x) = LmaxSim (x) ◦ TI p (x)→F (x)
(6)
maxSim
k=1 k
Dimensions (mm) 0.3516 0.3516 0.4688 0.9375 0.5804
Algorithm 1 Atlas Construction for Intra-Sequence Setting 1: Load set of atlases A = {a1 , a2 , ..., an } where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. 2: Calculate similarity between each pair of atlas images Sim(Ii (x), Ij (x)) = Corr(Ii (x), Ij (x)) 3: Using a K-means clustering function, cluster atlases into N classes Ac = {Ac1 , Ac2 , ..., AcN } where Ac = K means(M DS(Corr(Ii (x), Ij (x)))) p 4: For each class of atlases, create probabilistic atlas Ij (x) p and label Lj (x): m Pj k Ij (x) ◦ TIjk (x)→IjR (x) Ijp (x) = m1j
k=1
Lpj (x) =
512 × 512 × 140 512 × 512 × 70 512 × 512 × 80 256 × 256 × 100 448 × 448 × 50
In Equations 3 and 4, I and L are an image in the given class and its label, respectively, mj is the number of atlases in the class, and I p and Lp are the probabilistic atlas image and label, respectively. I R is a common reference image defined by minimizing a stack entropy cost function, which is the entropy of the stack of input images in each spatial location. The pixel value in each spatial location of the reference image is defined such that the cost function is minimized [19]. Algorithm 1 summarizes the atlas construction steps for the intra-sequence setting.
C. Intra-Sequence Setting: Results To evaluate the performance of the proposed algorithm, the segmentation results for the entire volume were compared to the ground-truth results (i.e., manual segmentation) using Dice Similarity Coefficient (DSC)6 and the point-to-surface distance error7 . A 5-fold cross-validation was performed with 50 and 6 For
2|A∩B|
two sets A and B, DSC is defined as |A|+|B| . calculate the point-to-surface distance error, first, for each point location in the result, a streamline is created in the normal direction by extracting both the point coordinates and the surface Normal. Next, the point-to-surface distance error is calculated by measuring the streamline distance between the point of origin in the resulting surface and the point at which the streamline intersects the ground-truth surface. 7 To
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
5
350 datasets, selected randomly, used for testing and atlas construction, respectively. The best configuration of the algorithm (i.e., 10 classes of atlases) yielded mean DSC of 93.16% ± 3.67%. The median value for the DSC was 93.98%. Figure 1 shows the distribution of the segmentation accuracy (DSC) for all cases. The mean point-to-surface distance error was 2.86 mm. Figure 2 shows the distribution of the point-to-surface distance errors for all cases. The processing time to segment a target 3D MR image was 2 min on a dual core Intel processor (Intel(R) Core(TM) i7-3770 CPU 3.40 GHz).
Fig. 3. Intra-sequence setting: accuracy (mean DSC) vs. number of atlas classes.
error of 0.07 mm for all test cases (2.93 mm to 2.86 mm)8 and 2.0 mm for the first decile of test cases (7 mm to 5 mm). The Wilcoxon signed-rank test was run on the DSC results of 1 class and 10 classes to confirm that they were statistically different (P < 0.001).
Fig. 1. Intra-sequence setting: DSC value comparison between the segmentation and ground-truth results for Dixon images
Fig. 4. Intra-sequence setting: accuracy (mean DSC) vs. number of atlas classes for first decile of test cases.
Fig. 2. Intra-sequence setting: point-to-surface distance errors (mm) between the segmentation and ground-truth results for Dixon images
One aspect of the proposed algorithm is the number of classes of atlases. Figure 3 shows the mean DSC values for the test images for different numbers of classes used to cluster the atlases. It is interesting to observe that there is an optimal number of classes of atlases (i.e., 10) that yields the best results (i.e., mean DSC of 93.16%). Although the increase in DSC for all test cases was small (1%: 92% to 93%), for the first 10% (decile) of test cases with lowest accuracy, the increase in accuracy was 3% (87% to 90%) (Figure 4); a decrease in point-to-surface distance
The effect of increasing the number of classes of atlases on the proposed algorithm performance can also be studied from the perspective of number of cases that meet a given threshold (e.g., 85%) for segmentation accuracy. Figure 5 shows the percentage of cases with DSC values below 85%. It is seen that increasing the number of classes of atlases from 1 to 10 decreases the percentage of failed cases from 7.6% to 2.8%; a decrease of 12 cases in our experiments of 250 test cases. We also investigated the effect of the number of atlas images used to create the probabilistic atlases on the accuracy of the proposed algorithm. We ran an experiment with one probabilistic atlas generated from different numbers of images 8 The 95% confidence intervals of point-to-surface distance error for one class of atlases and 10 classes of atlases for all cases were [2.84mm 3.02mm] and [2.79mm 2.93mm], respectively.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 6
IEEE TRANSACTIONS ON MEDICAL IMAGING
Fig. 5. Intra-sequence setting: percentage of cases with DSC below 85%.
(i.e., 1, 10, 20, 30, 40, and 50). It was found out that while increasing the number of images used to generate one probabilistic atlas from 1 to 10 increased the mean DSC significantly (83% to 90%), increasing the number of images further had insignificant effect on the accuracy. This means that in order to achieve high accuracy, a small number of atlas images is sufficient. In a previous work [30], the same dataset as the one used in this paper for the intra-sequence setting (section II-A) was segmented with a method based on a 3D local edge detection using phase congruency maps and Poisson surface reconstruction (Poisson Surface Reconstruction algorithm). The method also incorporated a mean shape information from a single atlas using a Laplacian framework. The reported mean DSC was 88% ± 4. The proposed algorithm in this paper outperforms the previous work in terms of accuracy of results (i.e., 93% vs. 88%). The main difference between the two algorithms are that, for the intra-sequence setting, our algorithm is a pure atlas-based segmentation method (with no edge information used from images) which combines atlas selection and probabilistic atlas approaches whereas the algorithm in [30] uses Poisson surface reconstruction to generate the 3D surface of the breast based on edge information coming from phase congruency maps. The algorithm in [30] uses a single atlas only to refine (or constrain) the 3D shape obtained from edgemaps. Figure 6 shows a sample test image, the corresponding bestmatch atlas, and the segmentation result for a single slice of Dixon images generated by the proposed algorithm for the intra-sequence setting. IV. P ROPOSED S EGMENTATION A LGORITHM FOR I NTER -S EQUENCE S ETTING In this section, the configuration of the proposed algorithm for the inter-sequence setting is presented where the atlas and target images are from different MR imaging sequences namely Dixon images and MR T1-weighted images, respectively. The main objective in the inter-sequence setting of the proposed algorithm is to create an atlas which is robust to intensity variation between different MR imaging sequences. Such features must represent information about the structures in the image. Although conventional edge detectors such as
Fig. 6. Intra-sequence setting (Dixon images). Top - Target image (left to right): original image, ground-truth result. Bottom - (left to right): best-match probabilistic atlas, segmentation result.
Canny [31] generate an edgemap of the image, they are not capable of generating features that are robust to variance in intensity or brightness. In contrast, phase congruency map (PCM) has been shown to be a robust image feature [32]. In a PCM, the local Fourier components of the image are all in phase (congruent) in locations where there are meaningful edges in the image. A PCM can be used to detect structural characteristics of an image in a way that is invariant to image intensity and robust to noise [32]. Equation 7 calculates the phase congruency of an image at location x (P C(x)) where Elocal (x) is the local energy of the image, T is a threshold to suppress the effect of noise on the local energy of the image at that location, Mn represents the amplitude of the nth Fourier component, and is set to a small number to avoid division by 0 [32]. The details of phase congruency implementation has been given in [30] where a monogenic signal is used. P C(x) =
bElocal (x) − T c Σn Mn (x) +
(7)
When the atlas and target images are from different MR imaging sequences, the intensity configuration of atlas and target images will be different potentially leading to poor registration results. To overcome this problem, we propose to use PCM to create atlas and perform atlas matching. This means that image features that are robust to intensity variations (e.g., real edges) are used to construct the atlas. Subsequently, when a target image arrives, its intensity invariant edges are computed, registered, and compared to that of the atlas images making it possible to segment a sequence of MR images different than that of the atlas images. A. Inter-Sequence Setting: Atlas Construction The atlas construction for the inter-sequence setting consists of clustering function and the probabilistic atlas building.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
1) Clustering Function: Similar to the intra-sequence setting (section III-A1), first, a (dis)similarity matrix of atlas images is created with the difference that in the inter-sequence setting, the cross-correlation measure is applied to PCMs of atlas images.
Sim(Ii (x), Ij (x)) = Corr(P CM (Ii (x)), P CM (Ij (x))) (8) A 2D distance map of PCMs of all images are generated using MDS algorithm and the K-means algorithm is used to cluster images into different classes based on the distance map of the corresponding PCMs. 2) Probabilistic Atlas: In order to create a probabilistic atlas for each class, first, the groupwise registration is applied to the PCMs of atlas images to generate the probabilistic PCM for each class (Equation 9). Next, the registration transformations resulted from groupwise registration of PCMs are applied to the labels in the class to generate the probabilistic label of the class (Equation 10). P CMjp (x)
mj 1 X P CM (Ijk (x)) ◦ TP CM (Ijk (x))→IjR (x) = mj k=1 (9)
7
p P CMmaxSim (x) = arg maxj f ine Corr(P CM (F (x)), P CMjp (x) ◦ TPafCM ) p (x)→P CM (F (x)) j
(11) Once the best-match class is found, the corresponding probabilistic PCM is registered (via nonrigid registration [29]) to the PCM of the target image. The registration transformation is then applied to the probabilistic label of the best-match class to create the segmentation result for the target image LA F (x). p nonrigid LA F (x) = LmaxSim (x) ◦ TP CM p
maxSim (x)→P CM (F (x))
C. Inter-Sequence Setting: Post-Processing Although our approach uses a robust atlas based on PCM, we anticipate that the registration result in the inter-sequence setting may require further refinement to achieve a higher accuracy. In general, to improve the accuracy of results, it is important to incorporate information taken directly from the target image to guide the atlas generated label. To do so, we use an active contour model [24] (ACM) in which the energy function E is minimized. Z1 Z1 E(v) =
mj
Lpj (x) =
1 X k Lj (x) ◦ TP CM (Ijk (x))→IjR (x) mj
(10)
where I R is a common reference used by groupwise registration. Algorithm 2 summarizes the atlas construction steps for the inter-sequence setting. Algorithm 2 Atlas Construction for Inter-Sequence Setting 1: Load set of atlases A = {a1 , a2 , ..., an } where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. 2: Calculate similarity between each pair of atlas images Sim(Ii (x), Ij (x)) = Corr(P CM (Ii (x)), P CM (Ij (x))) 3: Using a clustering function, cluster atlases into N classes Ac = K means(M DS(Corr(Ii (x), Ij (x)))) 4: For each class of atlases, create probabilistic PCM P CMjp (x) and probabilistic label Lpj (x): m Pj P CMjp (x) = m1j P CM (Ijk (x)) ◦ TP CM (Ijk (x))→IjR (x) Lpj (x)
=
1 mj
m Pj k=1
[Eint (v) + Eext (v)]ds 0
k=1
k=1
(12)
(13)
0
where v(s) is a surface, v(s) = (x(s), y(s)), and Eint and Eext are the internal and external energies of the image, respectively. The internal energy is defined as [24]: Eint (v) =
X 1 X (α |vs |2 + β |vss |2 ) 2 s s
(14)
where vs and vss are the first and second derivatives of the surface v, respectively, and α and β are weighting parameters to constrain the degree of elasticity and rigidity of the surface, respectively. The external energy represents the edge information of the image, which is defined as: Eext (v) = −|∇(Gσ (x, y, z) ∗ I(x, y, z))|2
(15)
where Gσ and I are the Gaussian function and the image, respectively, and ∇ is the gradient operator. For our experiments, we used the implementation in [33] where α = 0.5, β = 0.1 and Eext is calculated using Canny edge detector [31]. We used the label generated by atlas, LA F (x) from section IV-B, as the initial contour to the ACM.
Lkj (x) ◦ TP CM (Ijk (x))→IjR (x)
B. Inter-Sequence Setting: Atlas Matching When a target image F (x) arrives, first, its PCM (P CM (F (x))) is created, aligned (via affine registration [29]), and compared to the probabilistic PCMs of all classes to find p the best-match atlas PCM (P CMmaxSim (x)).
D. Inter-Sequence Setting: Results To evaluate the performance of the proposed algorithm, the segmentation results for the entire volume were compared to the ground-truth results (i.e., manual segmentation) using DSC and the point-to-surface distance error. A total of 17 cases were available from 9 different volunteer subjects using different machines (Table I). The best configuration of the algorithm (i.e., 10 classes of atlases) yielded mean DSC of 84.60% ± 6.5%. The median value for the DSC was 86.90%. Figure 7
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 8
IEEE TRANSACTIONS ON MEDICAL IMAGING
shows the distribution of the segmentation accuracy (DSC) for all cases. The mean point-to-surface distance error was 2.96 mm. Figure 8 shows the distribution of the point-to-surface distance errors for all cases. The processing time to segment a target 3D MR image was 2 min and 30 s on a dual core Intel processor (Intel(R) Core(TM) i7-3770 CPU 3.40 GHz).
Fig. 7. Inter-sequence setting: DSC value comparison between the segmentation and ground-truth results for MR T1-weighted images
not show a significant difference, the point-to-surface distance error confirmed the corresponding results were significantly different (P
1
Automated Segmentation of Breast in 3D MR Images Using a Robust Atlas Farzad Khalvati, Cristina Gallego Ortiz, Sharmila Balasingham, Anne L. Martel
Abstract—This paper presents a robust atlas-based segmentation (ABS) algorithm for segmentation of the breast boundary in 3D MR images. The proposed algorithm combines the wellknown methodologies of ABS namely probabilistic atlas and atlas selection approaches into a single framework where two configurations are realized. The algorithm uses phase congruency maps to create an atlas which is robust to intensity variations. This allows an atlas derived from images acquired with one MR imaging sequence to be used to segment images acquired with a different MR imaging sequence and eliminates the need for intensity-based registration. Images acquired using a Dixon sequence were used to create an atlas which was used to segment both Dixon images (intra-sequence) and T1-weighted images (inter-sequence). In both cases, highly accurate results were achieved with the median Dice Similarity Coefficient (DSC) values of 94% ± 4% and 87 ± 6.5%, respectively. Index Terms—Breast MRI, Image Segmentation, Atlas-Based Segmentation, Image Registration.
I. I NTRODUCTION
M
AGNETIC resonance imaging (MRI) of breast is an increasingly common approach for monitoring and detection of breast cancer mainly due to higher sensitivity and no ionizing radiation compared to conventional x-ray mammography. MR imaging is now also recommended for screening women who are known to be at a higher risk of breast cancer [1]. The segmentation of breast tissue in MR images has a number of important applications. The removal of unwanted structures such as the chest wall and the heart makes it possible to improve the quality of 3D visualisation. In order to use MRI to assess breast density, a significant risk factor and an important biomarker in determining the possibility of breast cancer [2], it is essential to determine the total volume of the breast. In Computer-Aided Diagnosis (CAD), breast segmentation allows the software to ignore tissue outside of the breast, reducing the false positive rate. Finally, in order to plan therapeutic interventions such as radiotherapy using MRI, accurate segmentations of the target organs will be required [3]. Atlas-based segmentation (ABS) is a well established and widely used technique for extracting contours from medical images of different organs (e.g., prostate) [4], [5], [6], [7], [8], [9]. In this method, processed images are stored in a Copyright (c) 2014 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to [email protected]. F. Khalvati is with the Sunnybrook Research Institute, Toronto, Ontario, Canada. e-mail: [email protected] C. Gallego Ortiz, S. Balasingham, and A. L. Martel are with the Sunnybrook Research Institute, Toronto, Ontario, Canada.
database called a set of atlases. Each atlas consists of an image and its gold standard (i.e., manual segmentation or label). The atlas images are first registered to the target image and then the resulting transformations are applied to the labels. Finally, the deformed labels are combined to achieve the final segmentation result. In general, there are three approaches to design an ABS algorithm: multi-atlas, atlas selection, and probabilistic atlas. The multi-atlas approach registers a limited number of atlas images to the target image to create multiple (deformed) labels, which are then combined using a fusion algorithm to generate the final result. Different fusion algorithms for multi-atlas approaches have been proposed in the literature. Majority voting is the simplest approach where each pixel is assigned to the label which occurs with the highest frequency. More complex methods attempt to weigh the segmentation results from each atlas based on their estimated accuracy [5]. Langerak et al. [9] proposed a multi-atlas segmentation method for prostate MR images based on an iterative label fusion approach. First, all atlas images are registered to the target image to obtain a set of deformed labels. The deformed labels are fused together using a weighted majority voting method and then, the overlap of each label is calculated against the fused label. Labels with low overlaps are discarded and a new fused label is recalculated. As another label fusion method, the STAPLE algorithm [10] uses an expectation-maximization approach in which the unknown true segmentation and the performance of individual atlas segmentations are estimated simultaneously. The atlas selection approach proposed in [5] uses both STAPLE and majority voting as the post-processing stage to fuse the generated labels. Another approach is to assume that the performance of each segmentation is related to the similarity between the registered atlas image and the target image. In [11], it was found that weighting labels using a local image similarity measure given by the mean squared difference between target and registered atlas image gave better results than majority voting or STAPLE. Several atlas selection methods have been proposed which try to find the subset of atlas images that best represent a specific target image in order to improve accuracy. The mutual information (MI) between the atlas and the target image after registration has been used as a figure of merit; for example in [12] the target images are subdivided into regions and for each region the atlas image with the best local MI is selected, and in [5] the subset of atlas images that has an MI greater than a pre-set threshold level is used. Atlas selection approaches that require all atlas images to be
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 2
IEEE TRANSACTIONS ON MEDICAL IMAGING
registered to the target image are computationally prohibitive and limit the number of atlases that can be used. To reduce the computational cost of atlas selection, Aljabar et al. [13] proposed to register all atlas images and the target image to an arbitrary reference image and thus, the best-match atlas is selected by comparing the registered target image with the registered atlas images. In other words, the atlas registration can be performed off-line during the atlas construction phase reducing the computational time. In a recent work, Langerak et al. [14] proposed a method that clusters the atlas images using the pairwise registration results to group similar atlases together. For each cluster, an exemplar atlas is selected and these are then registered to the target image. Once the best exemplar is found, all of the atlas images belonging to that cluster are registered to the target image and the resulting labels are combined using a simple majority voting. Although the process of clustering the atlases off-line is extremely computationally expensive, the number of registrations required for each target image is somewhat reduced. The probabilistic atlas approaches are usually based on either sequential models or generative models. In sequential models [4], [7], the probabilistic atlas is generated by aligning (and averaging) multiple images and their labels in a dataset using a registration algorithm. In generative models [15], [16], [17], [18], the tissue classification and image deformation are performed simultaneously using an objective function optimized by deforming the tissue probability maps. The main difference between the two classes is that in the former, the images are first registered and then combined to generate a probabilistic atlas. Next, the resulting probabilistic atlas is registered to the target image and consequently, the probabilistic label is deformed via the registration transformation to generate the result for the target image. In contrast, in the latter model, the target image is segmented (classified) directly as a result of deforming the tissue maps of the image via optimizing the objective function. There are different variations of sequential models of probabilistic atlas algorithms in the literature. For instance, Martin et al. [4] proposed to create a probabilistic atlas for prostate by registering images to a manually picked reference image. Next, a mean image is created by averaging the registered images. Each image is registered against the mean image to produce a deformed label of the image. The deformed labels are then averaged to generate a probability map of the labels. To segment a target image, the mean image is registered to it and the probability map of the labels is deformed using the registration transformation. Dowling et al. [7] also proposed a similar approach. To build a probabilistic atlas using a sequential model, a reference image(s) which has already been manually segmented and selected plays a crucial role. In practice, however, generating a reference image is not straight forward; even with careful atlas selection, the selected image may not represent the population well enough and hence, the registration may not produce acceptable results. Balci et al. [19] extended the groupwise registration algorithm proposed in [20] to be an alternative to pair-wise registration when designing a sequential probabilistic atlas for segmenting brain MR images.
Groupwise registration aligns a set of images to a common reference image by generating a set of transformations that map the reference image to each of the images in the group. This makes the sequential probabilistic atlas more robust and unbiased to a particular reference image. Ashburner and Friston [16] proposed a generative model of probabilistic atlas which unifies tissue classification, bias correction, and template registration for brain MR images. The objective function is based on a mixture of Gaussian (MOG) models where the tissue classes are represented by Gaussian probabilities and the bias, which is a smooth and spatially varying artifact that corrupts the image, is taken into account by extra parameters in MOG model. In this method, a template (spatial priors) is used in the process of optimizing the objective function where the deformation (registration) is updated during the process. Recently, Iglesias et al. [17] proposed a framework for multi-atlas segmentation of brain in MR images. The proposed algorithm is an extension of the work presented in [16] where in contrast to [16], it builds multiple probabilistic atlases. Although the ABS algorithms have been applied extensively to different organs such as brain and prostate, there have been limited work on breast segmentation using ABS in the literature. Gubern et al. [21] proposed a probabilistic atlas algorithm to segment different structures in breast MR (T1-weighted) images. In this method, the probabilistic atlas is created using pair-wise registration of each atlas with a reference image. The segmentation is done in a Bayesian framework where for a given target image voxel, the segmentation is estimated by searching for the most probable label [21]. The generative models of probabilistic ABS algorithms are usually computationally expensive since they integrate multiple optimizations [16]. The main drawback of the sequential models of probabilistic ABS algorithms is that they usually depend on the image intensity, rather than image key characteristics. As a result, if intensity variations exist between the target image and atlas image, which is common when the two images are from different scanners, the ABS algorithms may perform poorly resulting in inaccurate segmentation results. Moreover, in the conventional ABS algorithms, the atlas and target images must be from the same image sequence (e.g., MRI T1-weighted) in order to obtain reasonable segmentation results. One possible approach to make the sequential models of probabilistic ABS intensity invariant is to use a registration algorithm that is robust to intensity variations. Heinrich et al. [22] proposed an intensity-robust feature for image registration based on modality independent neighbourhood descriptor, which can be used for nonrigid registration across different modalities and different sequences of one modality (e.g., MRI). In this paper, we propose an improved ABS algorithm which is novel in that it combines the atlas selection method with a probabilistic atlas generated using a sequential model. The atlas images are grouped into classes based on their pairwise similarities. Each class of atlases is then used to create a probabilistic atlas using a groupwise registration. The proposed algorithm is presented in two configurations
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
for the intra-sequence and inter-sequence settings where in the former, the atlas and target images are from the same MR imaging sequence and in the latter, they are from different MR imaging sequences. In order to tackle the intensity-dependence deficiency of conventional ABS algorithms, for the intersequence setting, we use phase congruency maps, which are robust to intensity variations, to build the probabilistic atlas1 . These intensity-invariant maps are also used to register the atlas and target images. In order to refine the atlas results, the actual edge information from the target image is also incorporated directly into the segmentation pipeline by using an active contour model. To reduce the computational time, our proposed approach performs the atlas selection in pre-nonrigid registration stage and hence, to segment a target image, it only requires one nonrigid registration operation. II. M ATERIALS The proposed algorithm in this paper has been configured for two settings: intra-sequence and inter-sequence. In the intra-sequence setting, the atlas and test (target) images2 are from the same MR imaging protocol namely Dixon images (section II-A). In the inter-sequence setting, the atlas and test images are from different MR imaging protocols namely Dixon images and MR T1-weighted images, respectively (section II-B). In the following sections, the datasets used in this study are described in more detail. A. Intra-Sequence Materials For the intra-sequence setting, the atlas and test data, used from a previous study [23], contained 400 breast MRI datasets (94×94×44 pixels, wtih a resolution of 2.56×2.56×2.56 mm) acquired using a modified Dixon sequence. The population consisted of two age groups; 320 women aged 15-30 and 80 women aged 40-60 who were healthy and recruited between December 2003 and December 2007. The fat and water images of breast were obtained on a 1.5T Signa Cvi MR system (GE, Waukesha, WI) using a dedicated coil for MR breast (Medical Advances, Milwaukee, WI). The subjects were scanned in the prone position and the images were acquired in the sagittal plane with a slightly modified version of the GE FSE Dixon sequence. The images used in this paper were of zero phase shift. Out of 400 MRI datasets, 350 datasets were selected randomly to create the atlas and 50 datasets were set aside and used to evaluate the segmentation performance. A k-fold cross-validation was performed with k = 5 and the results were averaged. The images were contoured by an expert user by manual initialization of an active contour technique [24] where the segmentation results were corrected by the user to create ground-truth results. B. Inter-Sequence Materials For the inter-sequence setting, the data used to create the atlas was the same as the one used for the intra-sequence 1 In the remaining of the paper, probabilistic atlas refers to a sequential model of the probabilistic atlas approach. 2 In the remaining of the paper, test image and target image are used interchangeably.
3
setting as described in the previous section, which contained 350 breast MRI datasets (94 × 94 × 44 pixels) of Dixon imaging sequence. The test data comprised MR T1-weighted images acquired in 2013 using 3 different MRI scanners at the Sunnybrook Health Sciences Centre, Toronto, Ontario, Canada, for 9 healthy volunteer subjects. The age of the subjects ranged from 33 to 60 years and their BMIs ranged between 20.1 to 30.1 kg/m2 . All data was obtained under the local institutional research ethics board. For 8 subjects, both breasts were scanned while for one subject, only the right breast was scanned. The subjects were scanned in the prone position and the images were acquired in the sagittal plane. In order to evaluate the performance of the algorithm on different scanners, three scanners of both 1.5T and 3T magnets were used with different image dimensions. Table I summarizes the scanners and acquired data used to evaluate the performance of the proposed algorithm. The ground-truth results were created by manual segmentation of the images using TurtleSeg [25] (Interactive 3D Image Segmentation) software. III. P ROPOSED S EGMENTATION A LGORITHM FOR I NTRA -S EQUENCE S ETTING Our ABS algorithm combines a sequential model of probabilistic atlas and atlas selection approaches to offer an efficient ABS method. The goal is to increase the accuracy of the results by diversifying the atlas. In this section, a configuration of the algorithm for the intra-sequence setting is presented where the atlas and target images are from the same MR image sequence. A. Intra-Sequence Setting: Atlas Construction Atlas construction consists of two phases: clustering function and the probabilistic atlas building. 1) Clustering Function: Let F (x) be a target 3D image (or a Fixed image) and A = {a1 , a2 , ..., an } be a set of atlases where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. The first step is to cluster the atlases into N 3 separate classes so that similar atlases are grouped into one class. This requires a similarity measure that calculates the match between any given pair of atlas images in A, Sim(Ii (x), Ij (x)). This builds a (dis)similarity matrix of atlas images where the matrix represents pairwise dissimilarities between all atlas images. For the similarity measure, we use the cross-correlation coefficient4 . Sim(Ii (x), Ij (x)) = Corr(Ii (x), Ij (x))
(1)
For n atlas images, the similarity function Sim creates (dis)similarity values, which are fed to a multidimensional scaling (MDS) algorithm [26] to create a 2D distance map of all images5 . Next, the K-means clustering algorithm [27] is used to cluster images into N classes based n(n−1) 2
3 As it will be discussed in section III-C, for a given training dataset, the optimal value for N can be determined by running experiments across different values of N . 4 Mutual information may also be used as the similarity measure. 5 It is assumed that the atlas images are roughly aligned. Otherwise, an affine registration must be used to align the images before clustering.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 4
IEEE TRANSACTIONS ON MEDICAL IMAGING
TABLE I D ESCRIPTION OF THE BREAST MR T1- WEIGHTED IMAGES USED FOR INTER - SEQUENCE SETTING Total studies 8 4 2 2 1
Scanner GE MR7250 3T, GE MR7250 3T, Phillips Phillips
Dimensions (pixels)
GE Signa CV/i 1.5T, 8 channel coil Sentinelle 16 channel Vanguard Breast MRI tabletop coil Sentinelle 16 channel Vanguard Breast MRI tabletop coil Achieva 3T, Philips 7 channel Sense Breast coil Achieva 3T, Philips 7 channel Sense Breast coil
on the distance map where each class of atlases Acj has mj atlases where j ∈ [1, N ]. Acj = {a1 , a2 , ..., amj }
(2)
The union and intersection of atlas classes will give the set of all atlases A and an empty set, respectively. 2) Probabilistic Atlas: At this point, each class consists of several atlas images with the corresponding labels. The next step in atlas construction is to create a probabilistic atlas for each class of atlases. To create the probabilistic atlas, we use a groupwise registration technique which aligns a set of images to a common reference image. Balci et al. [19] proposed a framework for groupwise registration of images which uses the stack entropy cost function and a multi-resolution B-spline nonrigid deformation to generate the set of image transformations. The motivation behind using stack entropy as a cost function was that by aligning the images accurately, the intensity values of pixels in the corresponding locations in the stacked images should not vary significantly, which means that the stack entropy should be low [19]. The proposed approach uses a combination of global and local transformations where the former is an affine transformation and the latter is a nonrigid deformation based on B-splines [28]. In order to create the probabilistic atlas for each class, first, the groupwise registration is applied to atlas images in each class to generate the probabilistic image for the class (Equation 3). The transformations are applied to labels of each class to generate the probabilistic label of the class (Equation 4). Ijp (x) =
mj 1 X k Ij (x) ◦ TIjk (x)→IjR (x) mj
(3)
mj 1 X k Lj (x) ◦ TIjk (x)→IjR (x) mj
(4)
k
× × × × ×
0.3516 0.3516 0.4688 0.9375 0.5804
× × × × ×
2 2 2 2 2
Lpj (x)
=
1 mj
k=1 m Pj
k=1
Lkj (x) ◦ TIjk (x)→IjR (x)
B. Intra-Sequence Setting: Atlas Matching In the atlas matching phase, the target image F (x) is first aligned (via affine registration [29]) and compared to the probabilistic atlas images to find the best-match atlas image p ImaxSim (x). p f ine ImaxSim (x) = arg maxj Corr(F (x), Ijp (x) ◦ TIaf ) p j (x)→F (x) (5) Once the best-match probabilistic atlas image is found, it is registered (via nonrigid registration [29]) against the target image to obtain the registration transformation. The obtained transformation is applied to the best-match atlas label to generate the segmentation result LA F (x). p nonrigid LA F (x) = LmaxSim (x) ◦ TI p (x)→F (x)
(6)
maxSim
k=1 k
Dimensions (mm) 0.3516 0.3516 0.4688 0.9375 0.5804
Algorithm 1 Atlas Construction for Intra-Sequence Setting 1: Load set of atlases A = {a1 , a2 , ..., an } where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. 2: Calculate similarity between each pair of atlas images Sim(Ii (x), Ij (x)) = Corr(Ii (x), Ij (x)) 3: Using a K-means clustering function, cluster atlases into N classes Ac = {Ac1 , Ac2 , ..., AcN } where Ac = K means(M DS(Corr(Ii (x), Ij (x)))) p 4: For each class of atlases, create probabilistic atlas Ij (x) p and label Lj (x): m Pj k Ij (x) ◦ TIjk (x)→IjR (x) Ijp (x) = m1j
k=1
Lpj (x) =
512 × 512 × 140 512 × 512 × 70 512 × 512 × 80 256 × 256 × 100 448 × 448 × 50
In Equations 3 and 4, I and L are an image in the given class and its label, respectively, mj is the number of atlases in the class, and I p and Lp are the probabilistic atlas image and label, respectively. I R is a common reference image defined by minimizing a stack entropy cost function, which is the entropy of the stack of input images in each spatial location. The pixel value in each spatial location of the reference image is defined such that the cost function is minimized [19]. Algorithm 1 summarizes the atlas construction steps for the intra-sequence setting.
C. Intra-Sequence Setting: Results To evaluate the performance of the proposed algorithm, the segmentation results for the entire volume were compared to the ground-truth results (i.e., manual segmentation) using Dice Similarity Coefficient (DSC)6 and the point-to-surface distance error7 . A 5-fold cross-validation was performed with 50 and 6 For
2|A∩B|
two sets A and B, DSC is defined as |A|+|B| . calculate the point-to-surface distance error, first, for each point location in the result, a streamline is created in the normal direction by extracting both the point coordinates and the surface Normal. Next, the point-to-surface distance error is calculated by measuring the streamline distance between the point of origin in the resulting surface and the point at which the streamline intersects the ground-truth surface. 7 To
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
5
350 datasets, selected randomly, used for testing and atlas construction, respectively. The best configuration of the algorithm (i.e., 10 classes of atlases) yielded mean DSC of 93.16% ± 3.67%. The median value for the DSC was 93.98%. Figure 1 shows the distribution of the segmentation accuracy (DSC) for all cases. The mean point-to-surface distance error was 2.86 mm. Figure 2 shows the distribution of the point-to-surface distance errors for all cases. The processing time to segment a target 3D MR image was 2 min on a dual core Intel processor (Intel(R) Core(TM) i7-3770 CPU 3.40 GHz).
Fig. 3. Intra-sequence setting: accuracy (mean DSC) vs. number of atlas classes.
error of 0.07 mm for all test cases (2.93 mm to 2.86 mm)8 and 2.0 mm for the first decile of test cases (7 mm to 5 mm). The Wilcoxon signed-rank test was run on the DSC results of 1 class and 10 classes to confirm that they were statistically different (P < 0.001).
Fig. 1. Intra-sequence setting: DSC value comparison between the segmentation and ground-truth results for Dixon images
Fig. 4. Intra-sequence setting: accuracy (mean DSC) vs. number of atlas classes for first decile of test cases.
Fig. 2. Intra-sequence setting: point-to-surface distance errors (mm) between the segmentation and ground-truth results for Dixon images
One aspect of the proposed algorithm is the number of classes of atlases. Figure 3 shows the mean DSC values for the test images for different numbers of classes used to cluster the atlases. It is interesting to observe that there is an optimal number of classes of atlases (i.e., 10) that yields the best results (i.e., mean DSC of 93.16%). Although the increase in DSC for all test cases was small (1%: 92% to 93%), for the first 10% (decile) of test cases with lowest accuracy, the increase in accuracy was 3% (87% to 90%) (Figure 4); a decrease in point-to-surface distance
The effect of increasing the number of classes of atlases on the proposed algorithm performance can also be studied from the perspective of number of cases that meet a given threshold (e.g., 85%) for segmentation accuracy. Figure 5 shows the percentage of cases with DSC values below 85%. It is seen that increasing the number of classes of atlases from 1 to 10 decreases the percentage of failed cases from 7.6% to 2.8%; a decrease of 12 cases in our experiments of 250 test cases. We also investigated the effect of the number of atlas images used to create the probabilistic atlases on the accuracy of the proposed algorithm. We ran an experiment with one probabilistic atlas generated from different numbers of images 8 The 95% confidence intervals of point-to-surface distance error for one class of atlases and 10 classes of atlases for all cases were [2.84mm 3.02mm] and [2.79mm 2.93mm], respectively.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 6
IEEE TRANSACTIONS ON MEDICAL IMAGING
Fig. 5. Intra-sequence setting: percentage of cases with DSC below 85%.
(i.e., 1, 10, 20, 30, 40, and 50). It was found out that while increasing the number of images used to generate one probabilistic atlas from 1 to 10 increased the mean DSC significantly (83% to 90%), increasing the number of images further had insignificant effect on the accuracy. This means that in order to achieve high accuracy, a small number of atlas images is sufficient. In a previous work [30], the same dataset as the one used in this paper for the intra-sequence setting (section II-A) was segmented with a method based on a 3D local edge detection using phase congruency maps and Poisson surface reconstruction (Poisson Surface Reconstruction algorithm). The method also incorporated a mean shape information from a single atlas using a Laplacian framework. The reported mean DSC was 88% ± 4. The proposed algorithm in this paper outperforms the previous work in terms of accuracy of results (i.e., 93% vs. 88%). The main difference between the two algorithms are that, for the intra-sequence setting, our algorithm is a pure atlas-based segmentation method (with no edge information used from images) which combines atlas selection and probabilistic atlas approaches whereas the algorithm in [30] uses Poisson surface reconstruction to generate the 3D surface of the breast based on edge information coming from phase congruency maps. The algorithm in [30] uses a single atlas only to refine (or constrain) the 3D shape obtained from edgemaps. Figure 6 shows a sample test image, the corresponding bestmatch atlas, and the segmentation result for a single slice of Dixon images generated by the proposed algorithm for the intra-sequence setting. IV. P ROPOSED S EGMENTATION A LGORITHM FOR I NTER -S EQUENCE S ETTING In this section, the configuration of the proposed algorithm for the inter-sequence setting is presented where the atlas and target images are from different MR imaging sequences namely Dixon images and MR T1-weighted images, respectively. The main objective in the inter-sequence setting of the proposed algorithm is to create an atlas which is robust to intensity variation between different MR imaging sequences. Such features must represent information about the structures in the image. Although conventional edge detectors such as
Fig. 6. Intra-sequence setting (Dixon images). Top - Target image (left to right): original image, ground-truth result. Bottom - (left to right): best-match probabilistic atlas, segmentation result.
Canny [31] generate an edgemap of the image, they are not capable of generating features that are robust to variance in intensity or brightness. In contrast, phase congruency map (PCM) has been shown to be a robust image feature [32]. In a PCM, the local Fourier components of the image are all in phase (congruent) in locations where there are meaningful edges in the image. A PCM can be used to detect structural characteristics of an image in a way that is invariant to image intensity and robust to noise [32]. Equation 7 calculates the phase congruency of an image at location x (P C(x)) where Elocal (x) is the local energy of the image, T is a threshold to suppress the effect of noise on the local energy of the image at that location, Mn represents the amplitude of the nth Fourier component, and is set to a small number to avoid division by 0 [32]. The details of phase congruency implementation has been given in [30] where a monogenic signal is used. P C(x) =
bElocal (x) − T c Σn Mn (x) +
(7)
When the atlas and target images are from different MR imaging sequences, the intensity configuration of atlas and target images will be different potentially leading to poor registration results. To overcome this problem, we propose to use PCM to create atlas and perform atlas matching. This means that image features that are robust to intensity variations (e.g., real edges) are used to construct the atlas. Subsequently, when a target image arrives, its intensity invariant edges are computed, registered, and compared to that of the atlas images making it possible to segment a sequence of MR images different than that of the atlas images. A. Inter-Sequence Setting: Atlas Construction The atlas construction for the inter-sequence setting consists of clustering function and the probabilistic atlas building.
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging KHALVATI et al.: AUTOMATED SEGMENTATION OF BREAST IN 3D MR IMAGES USING A ROBUST ATLAS
1) Clustering Function: Similar to the intra-sequence setting (section III-A1), first, a (dis)similarity matrix of atlas images is created with the difference that in the inter-sequence setting, the cross-correlation measure is applied to PCMs of atlas images.
Sim(Ii (x), Ij (x)) = Corr(P CM (Ii (x)), P CM (Ij (x))) (8) A 2D distance map of PCMs of all images are generated using MDS algorithm and the K-means algorithm is used to cluster images into different classes based on the distance map of the corresponding PCMs. 2) Probabilistic Atlas: In order to create a probabilistic atlas for each class, first, the groupwise registration is applied to the PCMs of atlas images to generate the probabilistic PCM for each class (Equation 9). Next, the registration transformations resulted from groupwise registration of PCMs are applied to the labels in the class to generate the probabilistic label of the class (Equation 10). P CMjp (x)
mj 1 X P CM (Ijk (x)) ◦ TP CM (Ijk (x))→IjR (x) = mj k=1 (9)
7
p P CMmaxSim (x) = arg maxj f ine Corr(P CM (F (x)), P CMjp (x) ◦ TPafCM ) p (x)→P CM (F (x)) j
(11) Once the best-match class is found, the corresponding probabilistic PCM is registered (via nonrigid registration [29]) to the PCM of the target image. The registration transformation is then applied to the probabilistic label of the best-match class to create the segmentation result for the target image LA F (x). p nonrigid LA F (x) = LmaxSim (x) ◦ TP CM p
maxSim (x)→P CM (F (x))
C. Inter-Sequence Setting: Post-Processing Although our approach uses a robust atlas based on PCM, we anticipate that the registration result in the inter-sequence setting may require further refinement to achieve a higher accuracy. In general, to improve the accuracy of results, it is important to incorporate information taken directly from the target image to guide the atlas generated label. To do so, we use an active contour model [24] (ACM) in which the energy function E is minimized. Z1 Z1 E(v) =
mj
Lpj (x) =
1 X k Lj (x) ◦ TP CM (Ijk (x))→IjR (x) mj
(10)
where I R is a common reference used by groupwise registration. Algorithm 2 summarizes the atlas construction steps for the inter-sequence setting. Algorithm 2 Atlas Construction for Inter-Sequence Setting 1: Load set of atlases A = {a1 , a2 , ..., an } where ai = (Ii (x), Li (x)); Ii and Li are the 3D atlas image and its manual label, respectively. 2: Calculate similarity between each pair of atlas images Sim(Ii (x), Ij (x)) = Corr(P CM (Ii (x)), P CM (Ij (x))) 3: Using a clustering function, cluster atlases into N classes Ac = K means(M DS(Corr(Ii (x), Ij (x)))) 4: For each class of atlases, create probabilistic PCM P CMjp (x) and probabilistic label Lpj (x): m Pj P CMjp (x) = m1j P CM (Ijk (x)) ◦ TP CM (Ijk (x))→IjR (x) Lpj (x)
=
1 mj
m Pj k=1
[Eint (v) + Eext (v)]ds 0
k=1
k=1
(12)
(13)
0
where v(s) is a surface, v(s) = (x(s), y(s)), and Eint and Eext are the internal and external energies of the image, respectively. The internal energy is defined as [24]: Eint (v) =
X 1 X (α |vs |2 + β |vss |2 ) 2 s s
(14)
where vs and vss are the first and second derivatives of the surface v, respectively, and α and β are weighting parameters to constrain the degree of elasticity and rigidity of the surface, respectively. The external energy represents the edge information of the image, which is defined as: Eext (v) = −|∇(Gσ (x, y, z) ∗ I(x, y, z))|2
(15)
where Gσ and I are the Gaussian function and the image, respectively, and ∇ is the gradient operator. For our experiments, we used the implementation in [33] where α = 0.5, β = 0.1 and Eext is calculated using Canny edge detector [31]. We used the label generated by atlas, LA F (x) from section IV-B, as the initial contour to the ACM.
Lkj (x) ◦ TP CM (Ijk (x))→IjR (x)
B. Inter-Sequence Setting: Atlas Matching When a target image F (x) arrives, first, its PCM (P CM (F (x))) is created, aligned (via affine registration [29]), and compared to the probabilistic PCMs of all classes to find p the best-match atlas PCM (P CMmaxSim (x)).
D. Inter-Sequence Setting: Results To evaluate the performance of the proposed algorithm, the segmentation results for the entire volume were compared to the ground-truth results (i.e., manual segmentation) using DSC and the point-to-surface distance error. A total of 17 cases were available from 9 different volunteer subjects using different machines (Table I). The best configuration of the algorithm (i.e., 10 classes of atlases) yielded mean DSC of 84.60% ± 6.5%. The median value for the DSC was 86.90%. Figure 7
0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2347703, IEEE Transactions on Medical Imaging 8
IEEE TRANSACTIONS ON MEDICAL IMAGING
shows the distribution of the segmentation accuracy (DSC) for all cases. The mean point-to-surface distance error was 2.96 mm. Figure 8 shows the distribution of the point-to-surface distance errors for all cases. The processing time to segment a target 3D MR image was 2 min and 30 s on a dual core Intel processor (Intel(R) Core(TM) i7-3770 CPU 3.40 GHz).
Fig. 7. Inter-sequence setting: DSC value comparison between the segmentation and ground-truth results for MR T1-weighted images
not show a significant difference, the point-to-surface distance error confirmed the corresponding results were significantly different (P
