- Email: [email protected]
Anatomical atlas in the context of head and neck radiotherapy and its use to automatic segmentation
Biomedical Signal Processing and Control 7 (2012) 447–455
Contents lists available at SciVerse ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Anatomical atlas in the context of head and neck radiotherapy and its use to automatic segmentation Adriane Parraga a,∗ , Benoit Macq b , Mathieu De Craene c a
Universidade Estadual do Rio Grande do Sul, Dept. Engenharia de Sistemas Digitais, Rua Santa Maria, 2300, Guaíba, RS, Brazil Université Catholique de Louvain, Communication and Rem. Sen. laboratory, Place du Levant 2, Belgium Universitat Pompeu Fabra, Center for Computational Imaging and Simulation Technologies in Biomedicine (CISTIB), Information and Communication Technologies Dept., Pg Circumvallacio 8, Barcelona, Spain b c
a r t i c l e
i n f o
Article history: Received 18 July 2011 Received in revised form 6 December 2011 Accepted 12 December 2011 Available online 30 January 2012 Keywords: Anatomical atlas Registration Segmentation Head and neck anatomy
a b s t r a c t In this paper we present a methodology to form an anatomical atlas based on the analysis of dense deformation fields recovered by the Morphons non-rigid registration algorithm. The methodology is based on measuring the bending energy required to register the whole database to a reference, and the atlas is the one image in the database which yields the smallest bending energy when taken as reference. The suitability of our atlas is demonstrated in the context of head and neck radiotherapy through its application to a database with thirty-one computed tomography (CT) images of the head and neck region. In head and neck radiotherapy, CT is the most frequently used modality for the segmentation of organs at risk and clinical target volumes. One challenge brought by the use of CT images is the presence of important artifacts caused by dental implants. The presence of such artifacts hinders the use of intensity averages, thus severely limiting the application of most atlas building techniques described in the literature in this context. The results presented in the paper show that our bending energy model faithfully represents the shape variability of patients in the head and neck region; they also show its good performance in segmentation of volumes of interest in radiotherapy. Moreover, when compared to other atlases of similar performance in automatic segmentation, our atlas presents the desirable feature of not being blurred after intensity averaging. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The introduction of intensity modulated radiation therapy (IMRT) into clinical practice allows a better control of dose distribution over tumoral areas and the reduction of normal tissues jeopardizing. An efficient application of this technology in clinical practice raises the issues of adequacy and accuracy of selecting and delineating the clinical target volumes as well as surrounding organs at risk (OAR) to be preserved. Such delineation is typically performed by trained experts and is an extremely time consuming task. Manual segmentation also raises the risk of introducing intra- and inter-rater variabilities. Atlas based segmentation is a well known paradigm to assist the radiologist in the segmentation task [9]. In this paradigm, shape and intensity characteristics are encoded in the atlas, which is warped on the patient under investigation by a spatial transformation. Using this warping, volumes of interest (VOI) defined in the atlas can be projected onto the patients’ coordinate system.
∗ Corresponding author. E-mail address: [email protected] (A. Parraga). 1746-8094/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.bspc.2011.12.001
Significant efforts have been directed toward the development of templates for atlas based segmentation, mainly for the human brain using magnetic resonance (MR) as the main modality [25,29,17]. Some methods use a single subject anatomy, as in [5,6,22], which may be unfit to cope with the complexity of the image data and the variability of the structures under study. In these cases the atlas will be biased towards the specific anatomy of the reference subject. The arbitrary choice of the atlas as a single image representing the population introduces bias [13] in the atlas. In [11] and in [13,16] an average model is created to eliminate this bias from the chosen image. A widely used atlas of the human brain is the MNI atlas, that is the standard atlas template in statistical parameter mapping (SPM) [28]. This atlas was constructed using spatial normalization by linear registration with 9 degrees of freedom. Linear registration does not compensate for local shape differences in the brain, which induces blurring, not only in the averaged MR template but also in the tissue distribution maps obtained by averaging segmentations over all subjects. Atlases built from average images are not suited as an average shape template for non-rigid atlas-to-subject image segmentation [28]. Recently atlases based segmentation in the context of Head and Neck CT images (computed tomography) have been proposed to
448
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
segment the lymph nodes [8], and organs at risk for radiotherapy using a multi-atlas scheme [9,24]. Also automatic dental segmentation to guide dental care in the context of intensity-modulated radiotherapy has been proposed [19]. In the former work, the authors stated that the average-intensity atlas was blurred around teeth owing to artifacts, which limited the capabilities of accurate non-rigid registration, and consequently segmentation, around dental structures. A single image Atlas, like proposed in this paper, will not suffer from these limitations. What we propose here is the use of the deformation field as the main feature to build an atlas. Atlas construction based on deformation field analysis has been an important topic of research [29]. More closely related to our work is the one in [11], where it is proposed to build an average shape and average intensity atlas from a set of subjects. In that work, one image is chosen as an initial reference from a set of five magnetic resonance images of the brain, and an iterative procedure based on deformation fields is applied to construct an image that would represent the shape and intensity average of the population. They have shown convergence for any choice of the reference image to the final average anatomic template. However, only five images were used for template construction; such limited database may not be enough to assess local shape differences or to generalize the fast convergence into an unbiased template. In this paper, contrary to [11], we do not perform intensity average. This is mainly motivated by the fact that the main modality in the context of head and neck radiotherapy is computed tomography (CT), in which dental implants may cause metal artifacts in the volumes, making impracticable to average intensities as performed in most proposed methods [28]. Also, instead of performing an iterative construction of the shape average, we propose a criterion and a procedure to find, among the original images in a database, the one that is closest to this shape average and use this image as the atlas. Our results, obtained on a database of thirty one CT images, show that the performance of this “closest to average” atlas in automatic segmentation is as good as the performance of the average itself. Moreover, this atlas does not suffer from blurring introduced by averaging and interpolation procedures. The paper is organized as follows. In the next section we briefly describe the non-rigid registration method, called Morphons, used in this work. In Section 3, we describe the atlas building methodology, as well the implementation details. The numerical results are presented in Section 4, along with an illustrative qualitative evaluation. An analysis of these results and their implications is given in Section 5. Finally, Section 6 presents concluding remarks, also pointing out future work. 2. Image registration
are CT volumes, the data tend to be quite noisy, containing artifacts that cause intensity variations. These artifacts are illustrated in Fig. 1, which shows slices of an axial and a sagittal view of one particular subject of our database. 2.1. The Morphons method The Morphons method is a non-rigid registration method using an iterative deformation scheme where the displacements estimations are found from local phase difference [30,14]. To find local phase information, a set of quadrature filters, each one sensitive to structures in a certain direction, is applied to both the fixed image and the atlas. The information on local phase allows the detection of local structures, such as lines and edges [7]. In order to detect large structures a multiresolution scheme is used. This multiresolution scheme involves an iterative accumulation of the deformation field under the influence of certainty measures, which are associated with the displacement estimates found in each iteration. The formulation of the Morphons method is detailed in the sequel. Let x = (x, y, z) denote a point in the space, I(x) denote an image and Q(x) denote the impulse response of a quadrature filter. When the quadrature filter is applied to an image, the output is given by q(x) = Q (x) × I(x) = A(x)ejϕ(x)
Since the impulse response of a quadrature filters is complex, the output has a magnitude, A(x), and a phase value, ϕ(x). Applying the quadrature filter to a pair of images IF (x) and IM (x) – representing the fixed and moving images respectively – results in q1 (x) = Q (x) × IF (x) = A1 (x)ejϕ1 (x) q2 (x) = Q (x) × IM (x) = A2 (x)ejϕ2 (x)
(2)
Now, compute the product of q1 (x) and the complex conjugate of q2 (x): q1 (x)q∗2 (x) = A1 (x)A2 (x)ejϕ(x)
(3)
where the phase difference is the argument of this product: ϕ(x) = ϕ1 (x) − ϕ2 (x). The Morphons method determines the local displacement estimate dı to be performed in a certain filter direction ı as being proportional to the local phase difference of the filter responses in that direction. That is, at each voxel dı ∝ ϕ(x). At each voxel, a displacement estimate is found for each filter in the filter set, that is, a displacement field dı is obtained for each ˆ ı . These fields are combined into one displacement filter direction n field by solving a least squares problem: min
d
Registration is the process of finding a transformation that best matches two images according to a measure of similarity. One image is the reference, which remains fixed during the registration process, while the other is deformed in the geometric space of the reference. The reference image is also called fixed or target image and the transformed image is called the moving image. For the atlas building, it is important to use very accurate registration methods. The non-rigid algorithm used in this application is the Morphons algorithm. This method was chosen based on our previous work [23] and also in [12], in which we have shown with quantitative results the capability of this non-rigid multi-modality registration method to segment regions of interest (ROI) in head and neck CT images. Morphons extracts special features from the intensity images (output of quadrature phase filters) and works on these features instead of the original intensity. This makes the Morphons algorithm appropriate for multi-modality registration and noisy applications. Even though in our application all images
(1)
T
ˆ ı d − dı )]2 [cı (n
(4)
ı
ˆ ı is the direction of filter where d is the sought displacement field, n ı, and cı is a certainty measure associated to the i-th direction. This certainty measure is usually taken as the magnitude of Eq. (3), that is, cı = A1 (x)A2 (x) [14]. Iterative accumulation of a dense deformation field is done under the influence of the certainty measures: da =
ca da + ck (da + dk ) ca + ck
(5)
where da indicates the updated accumulated deformation field, da is the accumulated field from the previous iteration and dk is the displacement estimate derived in the current iteration k. ca and ck are certainty estimates associated with the accumulated deformation field and the displacement estimates, respectively. ck is the sum of certainty values cı in Eq. (4) for all filter directions. ca is
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
449
Fig. 1. Slices of an axial view and a sagittal view of one image of the database showing metal artifacts from dental implants.
found by an accumulation procedure by using the certainty values as certainties of themselves: ca =
ca2 + ck2 ca + ck
(6)
This iteration is applied at each level of the multiresolution scheme. From one level to another, the deformation field is upsampled and accumulated until the final level is reached [30]. 3. Anatomical atlas The aim of atlas construction is to have a unique image/volume that represents a population’s anatomy. The main purpose of the atlas in head and neck radiotherapy is to automatically segment the zones at risk and also the regions with high probability of tumoral propagation in the fat tissues [10]. The union of all volumes of interests, including some error margin, will be designated in the following as clinical target volumes (CTVs). Several methods have been proposed to create a brain atlas, but the suitability of such atlases in other anatomical regions has neither been demonstrated nor validated. In fact, in our application to head and neck we encountered difficulties that do not exist in brain MR images; these difficulties severely limit the application of most proposed methods in head and neck CT. Firstly it is prohibitive to have CT images/volumes of normal subjects due to the nature of the exam. Another major problem is the presence of artifacts, which invalidates the use of intensity averaging techniques. Recently, an atlas methodology for constructing head and neck atlases was investigated using eight patients [3], based on [11], using three different non-rigid registration algorithms. The drawback of their work is the degree of blurriness presented in the final atlas [2]. 3.1. Atlas model Our atlas building consists in finding a subject that best represents anatomically a population. In order to chose this subject, one needs an unambiguous numerical criterion which indicates that a certain image is the most representative in size/shape of human anatomy, this image being then chosen as the atlas. Non-rigid registration accounts for the differences of subjects due to size/shape variation, enabling the deformation of each of the subjects to be quantitatively and qualitatively compared. Accordingly, we will use deformation fields as the main feature to discriminate the best subject to be the atlas. The subject for which the average of the deformation field is closest to zero, after having all other subjects being registered into him, is the optimal shape
template or atlas. This is the principle of our method, assuming that all the images in the database are geometrically aligned. Let us now introduce formally our methodology. A deformation field is a function that represents a correspondence between two different subjects. Let D(x) = (Dx (x), Dy (x), Dz (x)) represent a deformation field for the whole image/volume, Dx (x), Dy (x) and Dz (x) being its components in each spatial direction. Let i and j be the indices of two subjects in a database of N images. The displacements from subject j to subject i in each direction x, y and z are described respectively by Dxij (x), Dyij (x) and Dzij (x). This triplet defines the deformation field that registers subject j into subject i and represents how much a subject has to be deformed to have the same size/shape of the other. After all subjects j have been registered in a certain subject i, the average of all the deformation fields for subject i will be found. The average displacements in each direction x, y and z are, respectively:
Dxi (x) =
Dyi (x) =
Dzi (x) =
1 N−1
1 N−1
1 N−1
N
Dxij (x)
(7)
Dyij (x)
(8)
Dzij (x)
(9)
j=1,j = / i N j=1,j = / i N j=1,j = / i
and Di (x) = (Dxi (x), Dyi (x), Dzi (x)) is the average deformation field for the subject i. The subject i whose average of the deformation field is closest to zero is the optimal shape template. To implement this rationale, we need to measure the closeness to zero in a concise manner, since Di (x) is a very dense information. We do that by using the norm of the magnitude, as defined in the following equations.
|Di (x)| =
Dxi (x)2 + Dyi (x)2 + Dzi (x)2
1 |Di (x)| m×n×p m
Di =
n
(10)
p
(11)
x=1 y=1 z=1
where m, n and p are the image dimensions. The quantity Di thus defined represents the mean bending energy of the deformation field for subject i. It represents, in a single number, how far the subject i stands from the average anatomy. Hence,
450
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455 Table 1 D values (in voxels) for each patient of the database.
Fig. 2. Scheme that illustrates the atlas building process. See text for details.
i
Subject
Di
i
Subject
Di
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 31
Subject1 Subject2 Subject3 Subject4 Subject5 Subject6 Subject7 Subject8 Subject9 Subject10 Subject11 Subject12 Subject13 Subject14 Subject15 Subject31
1.79 3.59 1.07 2.12 2.62 2.09 2.28 4.45 2.23 1.41 2.12 2.16 1.42 2.94 1.32 5.31
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Subject16 Subject17 Subject18 Subject19 Subject20 Subject21 Subject22 Subject23 Subject24 Subject25 Subject26 Subject27 Subject28 Subject29 Subject30
4.32 2.36 2.51 1.61 1.96 1.72 1.48 4.28 1.87 1.80 1.65 2.24 1.62 1.49 1.50
the subject that requires the minimum bending energy is the optimal anatomical atlas A, that is: A = Ik ,
where
k = arg(minDi )
(12)
4.1. Database
i
The block diagram in Fig. 2 illustrates the atlas building process. In this scheme, let each gray box be a subject in a database of 6 images. The boxes are disposed by their size/shape in relation to the central image, this latter being the most central anatomy in this database. Since we do not know a priori which image has the most representative size/shape, each and everyone at a time will be the fixed image and all the 5 subjects left will be non-rigidly registered into him. Based on the left diagram of Fig. 2, subject S1 starts being the atlas candidate. The deformation field resulting from registering image S2 into S1 is denoted by D12 , registering image S3 into S1 yields the deformation field D13 , and so on. After all images have been registered to S1, Eqs. (7), (8) and (9) are calculated for the set of deformation fields D1j from all registrations pairs 1j. Then the Eq. (11) is computed for subject S1. Next, all steps described above have to be repeated considering S2 as the fixed image – or the atlas candidate – as shown at the right in Fig. 2. The procedure is repeated, each subject at a time being the atlas candidate. The subject i that yields the minimum value, as specified in Eq. (12), is the best shape template for that population and is the proposed atlas A. In order to illustrate our atlas construction procedure, we took a real image of a hand and artificially created two others to have an enlarged and a reduced version of it. Fig. 3 shows the real hand on top, the enlarged version in the left bottom and the reduced one in the right bottom. In this experience, after registering both modified hands at the original one, we calculated the mean bending energy of the deformation field, finding a value of D1 = 0.44. Registering the enlarged hand and the original one into the reduced hand, the mean bending energy provided a value of D2 = 8.42. Then the reduced hand and the original one were registered into the enlarged hand, providing a mean bending energy value of D3 = 6.48. The mean bending energy of the average deformation field was around 15–20 times smaller considering the original hand as the central anatomy. This corroborates the expected behavior of our method, showing that our metric is capable of finding the image with a central anatomy in the population.
4. Results We describe in the sequel the application of the proposed atlas building procedure to a real database of head and neck images, as well as the atlas’ use in automatic segmentation of OARs and CTVs.
A database of 31 patients of 3D CT images previously segmented by a radiologist has been used to build the atlas. The size of the CT volumes is of 256 × 256 × 128 voxels with a voxel size of 0.9765 mm × 0.9765 mm × 2.1093 mm. All images in this database are male patients with age between 50 and 70 years old.
4.2. Implementation Before performing Morphons registration, all images in the database have to be at same geometric space; so we firstly apply a rigid registration. The geometric alignment is obtained by minimizing the mutual information between the fixed and the moving image [20] using a simultaneous perturbation stochastic approximation (SPSA) [4]. The rigid transformation is a superposition of a 3-D rotation and a 3-D translation and the registration parameter is a six-component vector consisting of three rotation angles and three translation distances [18]. The rigid registration is done only once, putting all database at the same geometric space of any chosen subject. The choice of this subject can be arbitrary, because linear transformation such as the rigid registration preserves the anatomy of the moving image and consequently does not affect the result. The rigid registration algorithm was performed using the ITK environment (http://www.itk.org) and Morphons was performed using Matlab. A total of 31 rigid registrations and 31×30 non-rigid registrations has been performed. For each image i being the fixed in the registration process, we have calculated D from Eq. (11). Table 1 summarizes the numerical results – the D values for all 31 patients, which are given in unities of voxels. It is seen that subject 3 was found to give the minimal value of D = 1.07. This means that subject 3 requires the minimum mean bending energy for all database to be deformed into him, and that the average displacement for each voxel is a distance equivalent to 1.07 voxels. Fig. 4 shows an axial and sagittal view of subject 3. To illustrate that subject 3 has indeed a central anatomy regarding the size/shape in relation to others in the database, we took two subjects with high values of D: patient 31 with D31 = 5.31 and patient 16 with D16 = 4.32. These high values mean that these subjects have extreme anatomies, that is, anatomies that are significantly distant from the population’s average. This is indeed the case, as we can observe in the Fig. 5, which presents an axial and a sagittal view of patient 31, and in Fig. 6 which presents an axial and a sagittal view of patient 16. This means that in average, those subjects require higher mean bending energy.
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
Fig. 3. Database of hands. Top: original hand; left bottom: artificially enlarged hand; right bottom: artificially reduced hand.
Fig. 4. Slice of patient 3. Left: axial view; right: sagittal view.
Fig. 5. Slice of patient 31. Left: axial view; right: sagittal view.
451
452
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
Fig. 6. Slice of patient 16. Left: axial view; right: sagittal view.
4.3. Atlas based segmentation Next, we apply our atlas to the task of automatic segmentation of the regions of interest in head and neck radiotherapy and compare the results to the manual segmentation. The regions of interest were manually contoured in all images by a physician, including the right and left parotid glands, the spinal cord, the clinical target volume and also the body contour. From each contour a binary image (mask) has been created, setting ones inside the contour (ROI) and zeroing the remainder of the grayscale image. We use this manual segmentation performed by an expert for comparison purposes. The atlas based segmentation requires that the atlas itself has been previously segmented, so that binary images are available for each region of interest. The segmentation of a given subject’s image is then achieved in two steps. First the atlas is registered into the subject and then the resulting deformation field is applied to the binary masks of the atlas. The result of this process is a set of binary masks for the patient, each one representing one region of interest. To evaluate the segmentation done by the atlas, the deformed masks are compared to the manual segmentation done by the physician. The metric used to measure spatial intersection of two masks was the similarity index (SI) [31]: SI = 2 ·
|MP ∩ MA | |MP | + |MA |
(13)
where | · | is the number of non-zero voxels of the masks and SI ∈ [0, 1]. The SI was calculated between the masks originated from the manual segmentation – MP – and the masks originated from the automatic segmentation using our atlas – MA – as in Eq. (13). Table 2 presents the resulting SIs for the ROIs of selected subjects from the database.1 The procedure in [11] constructs iteratively an average shape of the population, whereas we only pick, among the available (real) images, the one that is closest to this average. To compare the atlas obtained with our strategy with the one described in their work, we unbiased our atlas using their method. This is achieved by applying the inverse average deformation field into the first image chosen as the atlas, creating a new atlas. Then all images in the database are registered into this new version of the atlas, and then again the inverse average deformation field is applied to the present atlas, creating again a new version of the atlas. This process is done iteratively until convergence is reached, that is, until the new atlas does not change significantly from the previous iteration. This is also similar to the method proposed in [15]. This scheme is illustrated in Fig. 7. In this figure, the symbol A indicates the atlas A proposed here and the boxes represent the
1 The non-filled lacunas in Table 2 correspond to the contours that were not available in the Dicom-RT original images.
other images in the database. At the first iteration, all subjects, denoted by S1 . . . S6 in Fig. 7, are registered into the atlas A. Then, the average deformation field is calculated and its inverse is applied to the atlas A, creating a new version of an atlas, denoted by G1 . These steps are repeated, but now registering all images from the database in G1 , creating a new estimation of the average (Guimond’s) atlas, G2 . The convergence of our atlas A to the Guimond’s atlas occurred after 6 iterations and we called it atlas G. To illustrate it qualitatively, an axial and a sagittal views are shown of Guimond’s atlas in Fig. 8. To help the visualization of the atlas G, the A is also presented in Fig. 9 with a blue contour placed at the bones, and this contour is overlayed in G to facilitate the visualization of anatomic differences between those images. After atlas G was created, it was used to segment the ROIs previously described. The results are summarized in Table 3, where the same SI metric was used. It can be seen in this table that the performance of both atlases in the automatic segmentation task is quite similar. Indeed, performing an analysis of variance (ANOVA) to these data with a level of significance of 5% established that the results obtained with the two different atlases present no statistical differences (p-values p 0.05 were found). 5. Discussion In reference [11] an average atlas for the brain is arrived at after three iterations. However they have used a small database of only five patients for building the average shape. A similar work has reported more iterations for a larger database [13]. In this paper, we achieved convergence to the average atlas G after six iterations, with the major differences that our database consists of thirty one subjects and the first image used to build the Guimond’s atlas was not a random image from our database. The size of the database regarding the population that the atlas is supposed to represent is an important issue. In this work we have used a reasonably large database (thirty one subjects) to represent a somewhat restricted population, namely only male adults, so we expect the database to be sufficiently representative. The mean bending energy of our atlas A is 1.07, meaning that each voxel must be displaced an average of 1.07 voxels when registering all the subjects of the database in the atlas. Applying an iterative procedure to search for the average anatomy results in an alternative atlas whose final mean bending energy, after six iterations, is 0.20. Although this resulting atlas is certainly close to the central anatomy of the population, this difference does not seem to be relevant. Indeed, it has been seen that the performance of the two atlases in automatic segmentation is not distinguishable, and that anatomical differences between them are hardly perceivable for the naked eye. Hence our results imply that, in the quest for the average anatomy, there is a limiar beyond which it is useless
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
453
Table 2 Similarity index for segmenting the regions of interest from head and neck images using the proposed atlas A versus manual delineation done by a physician. ROI
Pat4
Pat10
Pat13
Pat16
Pat17
Pat23
Pat29
Body Spinal cord Right parotid Left parotid Right CTV Left CTV
0.98 0.69 0.70 0.68 – 0.73
0.96 0.74 0.69 0.57 – 0.47
0.97 0.65 0.54 0.30 0.60 0.57
0.97 0.55 0.70 – 0.70 0.64
0.96 0.51 0.54 0.49 0.55 0.59
0.98 0.77 0.63 0.65 – 0.70
0.98 0.73 0.60 0.57 0.65 –
Fig. 7. Scheme illustrating 3 iterations of Guimond’s method in the proposed Atlas A, generating an unbiased Atlas G.
Fig. 8. Axial and sagittal views of the atlas Guimond at the sixth (final) iteration [11]. The same blue contour from Fig. 9 was overlayed this image to help the visualization of anatomic difference between the two images. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Fig. 9. Axial and sagittal views of the atlas A, also used as the reference image for constructing Guimond [11]. A blue contour was placed in the bone structure to help visualization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Table 3 Similarity index for segmenting the regions of interest of head and neck images using atlas G versus manual delineation done by a physician. ROI
Pat4
Pat10
Pat13
Pat16
Pat17
Pat23
Pat29
Body Spinal cord Right parotid Left parotid Right CTV Left CTV
0.97 0.68 0.66 0.60 – 0.70
0.96 0.74 0.68 0.54 – 0.47
0.97 0.64 0.54 0.31 0.56 0.56
0.96 0.54 0.68 – 0.64 0.54
0.96 0.52 0.56 0.48 0.57 0.57
0.97 0.73 0.55 0.63 – 0.66
0.97 0.70 0.60 0.60 0.68 –
454
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
to proceed. On the other hand, most methods of atlas building use some sort of averaging and/or interpolation, which induces blurring; it is thus advisable not to proceed any further than necessary in such methods. In our database this limiar is slightly above 1 voxel and any image whose mean bending energy is smaller than this is very likely to be a good atlas. Whether this particular value can be taken as a general rule is a topic of further investigation. Recent results strongly indicate that the strategy of multi-atlas based segmentation is more effective than a single atlas to cope with the anatomy variability [9,24]. In the former work, clusters are created and then an average image of each cluster is an atlas. The methodology for atlas building proposed here could be used to create the cluster based on the mean bending energy for a multiatlas strategy. Instead of creating an average image, each cluster would have a median image chosen by our criterion that would be neither biased nor blurred. Also, a new image of a patient to be segmented would select the closest atlas based on the least mean bending energy. Our atlas could also be used as the initialization for a more sophisticated method of segmentation such as the active contour based segmentation as in [8], in which they have chosen an arbitrary image to be the atlas. 6. Conclusions and final remarks In this paper we have proposed a methodology to build an atlas that represents the size and shape variability present in human anatomy based on dense deformation fields using Morphons. In our previous work we have shown that Morphons is an effective strategy for performing non-rigid registration of the head and neck anatomy. Patients with high mean bending energy – a large D value – have more difference in anatomy in relation to the others, as we have shown. The atlas proposed here consists of the subject with the lowest D value within a population. The main advantage of our methodology is that the chosen atlas is an original image/volume and thus does not have any degree of blurriness imposed by interpolation or iterative processes, which could impose problems in future registration when using it for atlas based segmentation. On the other hand, it may happen, particularly when a restricted database is considered, that even the lowest D is still not small enough – meaning that the subject closest to the average anatomy is not close enough to it in order to be a good atlas. In such cases our atlas could be used as a complement to iterative procedures, by providing the best initial condition for them. One of our future work objectives is to improve the SI for some regions of interest being investigated. The goal is to build a statistical atlas based on [27] and also investigate the use of segmentation methods based on active contours [8], which the atlas based segmentations will be used as a initialization. Alternative metrics, such as in [1], will also be considered in future works. Finally we need to mention the computational complexity of our method. The number of registrations to be computed is N · (N − 1). Although this method is computationally intensive this does not pose any serious problems, since our calculations are performed off-line and all the registrations has to be done only once. Acknowledgements We gratefully acknowledge the support from the SIMILAR Network of Excellence, the Brazilian Education Ministry (MEC-CAPES) and the National Council of Technological and Scientific Development (CNPq) for financial support. The authors are grateful to Dr. Vincent Grégoire, Xavier Geets and Pierre Castadot for providing the images and their manual delineations. The authors wish to thank Vincent Nicolas for developing efficient visualization tools (http://www.medicalstudio.org/).
References [1] V. Arsigny, O. Commowick, X. Pennec, N. Ayache, A log-Euclidean framework for statistics on diffeomorphisms, in: Proceedings of the 9th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI’06), 2–4 October, Part I. No. 4190, 2006. [2] O. Commowick, V. Grégoire, G. Malandain, Atlas-based delineation of lymph node levels in head and neck computed tomography images, Radiother. Oncol. 87 (2) (2008) 281–289. [3] O. Commowick, G. Malandain, Evaluation of atlas construction strategies in the context of radiotherapy planning., in: Proceedings of the SA2PM Workshop (From Statistical Atlases to Personalized Models), Copenhagen, held in conjunction with MICCAI 2006, October 2006. [4] M.D. Craene, Dense deformation field estimation for pairwise and multisubjects registration, Ph.D. thesis, Université Catholique de Louvain, B - 1348 Louvain-la-Neuve Belgique, August 2005. [5] M.B. Cuadra, C. Pollo, A. Bardera, O. Cuisenaire, J.-G. Villemure, J.-P. Thiran, Atlas-based segmentation of pathological MR brain images using a model of lesion growth, IEEE Trans. Med. Imaging 23 (10) (2004) 1301–1314. [6] B.M. Dawant, S.L. Hartmann, J.-P. Thirion, F. Maes, D. Vandermeulen, P. Demaerel, Automatic 3D segmentation of internal structures on the head in MR images using a combination of similarity and free form transormations: Part I, metholody and validation on normal subjects, IEEE Trans. Med. Imaging 18 (10) (1999) 909–916. [7] W.T. Freeman, E.H. Adelson, The design and use of steerable filters, IEEE Trans. Pattern Anal. Mach. Intell. 13 (9) (1991) 891–906. [8] S. Gorthi, V. Duay, N. Houhou, M. Bach Cuadra, U. Schick, M. Becker, A.S. Allal, J.P. Thiran, Segmentation of head and neck lymph node regions for radiotherapy planning, using active contour based atlas registration, IEEE J. Sel. Top. Signal Process. 3 (2009) 135–147. [9] S. Gorthi, M. Bach Cuadra, U. Schick, P.-A. Tercier, A.S. Allal, J.-P. Thiran, Active multi-atlas based segmentation of head and neck CT images using active contour framework, in: MICCAI workshop on 3D Segmentation Challenge for Clinical Applications, 2011, p. 2010. [10] V. Grégoire, P. Levendag, J. Bernier, M. Braaksma, V. Budach, C. Chao, E. Coche, J. Cooper, G. Cosnard, A. Eisbruch, S. El-sayed, B. Emami, C. Grau, M. Hamoir, N. Lee, P. Maingon, K. Muller, H. Reychler, CT-based delineation of lymph node levels and related CTVs in the node-negative neck: DAHANCA, EORTC, GORTEC, NCIC,RTOG consensus guidelines, Radiother. Oncol. 69 (3) (2003) 227–236. [11] A. Guimond, J. Meunier, J.-P. Thirion, Average brain models: a convergence study, Comput. Vision Image Understanding 77 (2) (2000) 192–210. [12] G. Janssens, J. Orban de Xivry, B. Macq, Evaluation of nonrigid registration models for interfraction dose accumulation in radiotherapy, Med. Phys. 36 (9) (2009) 4268–4276. [13] S. Joshi, B. Davis, M. Jomier, G. Gerig, Unbiased diffeomorphic atlas construction for computational anatomy, Neuroimage 23 (Suppl. 1) (2009). [14] H. Knutsson, M. Andersson, Morphons: segmentation using elastic canvas and paint on priors, in: IEEE International Conference on Image Processing (Proceedings of the IEEE International Conference on Image Processing’05), Genova, Italy, September 2005. [15] P. Kochunov, J.L. Lancaster, P. Thompson, R. Woods, J. Mazziotta, J. Hardies, P. Fox, Regional spatial normalization: toward an optimal target, J. Comput. Assist. Tomogr. 25 (5) (2001) 805–816. [16] P. Lorenzen, M. Prastawa, B. Davis, G. Gerig, E. Bullitt, S. Joshi, Multi-modal image set registration and atlas formation, Med. Image Anal. 10 (June (3)) (2006) 440–451. [17] M.D. Craene, M.D. Duay, V. Macq, B. Pollo, C.J. Thiran, Dense deformation field estimation for atlas-based segmentation of pathological MR brain images, Comput. Methods Programs Biomed. (August) (2006) 66–75. [18] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Trans. Med. Imaging 16 (2) (1997) 187–198. [19] P. Maingon, L. Ramus, G. Odin, V. Grégoire, V. Darcourt, P.Y. Marcy, N. Guevara, M.H. Orlanducci, S. Marcie, G. Poissonnet, A. Bozec, O. Dassonville, F. Demard, L. Castillo, J. Santini, G. Malandain, DENTALMAPS: Automatic Dental Delineation for Radiotherapy Planning in Head and Neck Cancer, Int. J. Radiat. Oncol. Biol. Phys. (2011). [20] D. Mattes, D.R. Haynor, H. Vesselle, T.K. Lewellen, W. Eubank, Pet-ct image registration in the chest using free-form deformations, IEEE Trans. Med. Imaging 22 (1) (2003) 120–128. [22] H. Park, P. Bland, C.R. Meyer, Construction of an abdominal probabilistic atlas and its application in segmentation, IEEE Trans. Med. Imaging 22 (4) (2003) 483–492. [23] A. Parraga, A. Susin, J. Pettersson, B. Macq, M.D. Craene, Quality assessment of non-rigid registration methods for atlas-based segmentation in head-neck radiotherapy, in: Proceedings of IEEE International Conference on Acoustics, Speech, & Signal Processing, Honolulu, Hawaii, USA, April 2007. [24] L. Ramus, Conception et utilisation d’atlas anatomiques pour la segmentation automatique: application à la radiothérapie des cancers ORL, Ph.D. thesis, Université de Nice Sophia Antipolis, July 2011. [25] A. Rao, R. Chandrashekara, G.I. Sanchez-Ortiz, R. Mohiaddin, P. Aljabar, J.V. Hajnal, B.K. Puri, D. Rueckert, Spatial transformation of motion and deformation fields using nonrigid registration, IEEE Trans. Med. Imaging 23 (9) (2004) 1065–1076.
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455 [27] D. Rueckert, A.F. Frangi, J.A. Schnabel, Automatic construction of 3D statistical deformation models of the brain using non-rigid registration, IEEE Trans. Med. Imaging 22 (8) (2003) 1014–1025. [28] D. Seghers, E. D’Agostino, F. Maes, D. Vandermeulen, P. Suetens, Construction of a brain template from mr images using state-of-the-art registration and segmentation techniques., in: MICCAI (1), 2004, pp. 696–703. [29] Q. Wang, D. Seghers, E. D’Agostino, F. Maes, D. Vandermeulen, P. Suetens, A. Hammers, Construction and validation of mean shape atlas templates for
455
atlas-based brain image segmentation., in: Proceedings Information Processing in Medical Imaging – IPMI, vol. 3565, 2005, pp. 689–700. [30] A. Wrangsjö, J. Pettersson, H. Knutsson, Non-rigid registration using morphons., in: Proceedings of the 14th Scandinavian Conference on Image Analysis (SCIA’05), Joensuu, 2005 June. [31] A. Zijdenbos, B. Dawant, R. Margolin, A. Palmer, Morphometric analysis of white matter lesions in MR images: method and validation, IEEE Trans. Med. Imaging 13 (December (4)) (1994) 716–724.
Contents lists available at SciVerse ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Anatomical atlas in the context of head and neck radiotherapy and its use to automatic segmentation Adriane Parraga a,∗ , Benoit Macq b , Mathieu De Craene c a
Universidade Estadual do Rio Grande do Sul, Dept. Engenharia de Sistemas Digitais, Rua Santa Maria, 2300, Guaíba, RS, Brazil Université Catholique de Louvain, Communication and Rem. Sen. laboratory, Place du Levant 2, Belgium Universitat Pompeu Fabra, Center for Computational Imaging and Simulation Technologies in Biomedicine (CISTIB), Information and Communication Technologies Dept., Pg Circumvallacio 8, Barcelona, Spain b c
a r t i c l e
i n f o
Article history: Received 18 July 2011 Received in revised form 6 December 2011 Accepted 12 December 2011 Available online 30 January 2012 Keywords: Anatomical atlas Registration Segmentation Head and neck anatomy
a b s t r a c t In this paper we present a methodology to form an anatomical atlas based on the analysis of dense deformation fields recovered by the Morphons non-rigid registration algorithm. The methodology is based on measuring the bending energy required to register the whole database to a reference, and the atlas is the one image in the database which yields the smallest bending energy when taken as reference. The suitability of our atlas is demonstrated in the context of head and neck radiotherapy through its application to a database with thirty-one computed tomography (CT) images of the head and neck region. In head and neck radiotherapy, CT is the most frequently used modality for the segmentation of organs at risk and clinical target volumes. One challenge brought by the use of CT images is the presence of important artifacts caused by dental implants. The presence of such artifacts hinders the use of intensity averages, thus severely limiting the application of most atlas building techniques described in the literature in this context. The results presented in the paper show that our bending energy model faithfully represents the shape variability of patients in the head and neck region; they also show its good performance in segmentation of volumes of interest in radiotherapy. Moreover, when compared to other atlases of similar performance in automatic segmentation, our atlas presents the desirable feature of not being blurred after intensity averaging. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The introduction of intensity modulated radiation therapy (IMRT) into clinical practice allows a better control of dose distribution over tumoral areas and the reduction of normal tissues jeopardizing. An efficient application of this technology in clinical practice raises the issues of adequacy and accuracy of selecting and delineating the clinical target volumes as well as surrounding organs at risk (OAR) to be preserved. Such delineation is typically performed by trained experts and is an extremely time consuming task. Manual segmentation also raises the risk of introducing intra- and inter-rater variabilities. Atlas based segmentation is a well known paradigm to assist the radiologist in the segmentation task [9]. In this paradigm, shape and intensity characteristics are encoded in the atlas, which is warped on the patient under investigation by a spatial transformation. Using this warping, volumes of interest (VOI) defined in the atlas can be projected onto the patients’ coordinate system.
∗ Corresponding author. E-mail address: [email protected] (A. Parraga). 1746-8094/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.bspc.2011.12.001
Significant efforts have been directed toward the development of templates for atlas based segmentation, mainly for the human brain using magnetic resonance (MR) as the main modality [25,29,17]. Some methods use a single subject anatomy, as in [5,6,22], which may be unfit to cope with the complexity of the image data and the variability of the structures under study. In these cases the atlas will be biased towards the specific anatomy of the reference subject. The arbitrary choice of the atlas as a single image representing the population introduces bias [13] in the atlas. In [11] and in [13,16] an average model is created to eliminate this bias from the chosen image. A widely used atlas of the human brain is the MNI atlas, that is the standard atlas template in statistical parameter mapping (SPM) [28]. This atlas was constructed using spatial normalization by linear registration with 9 degrees of freedom. Linear registration does not compensate for local shape differences in the brain, which induces blurring, not only in the averaged MR template but also in the tissue distribution maps obtained by averaging segmentations over all subjects. Atlases built from average images are not suited as an average shape template for non-rigid atlas-to-subject image segmentation [28]. Recently atlases based segmentation in the context of Head and Neck CT images (computed tomography) have been proposed to
448
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
segment the lymph nodes [8], and organs at risk for radiotherapy using a multi-atlas scheme [9,24]. Also automatic dental segmentation to guide dental care in the context of intensity-modulated radiotherapy has been proposed [19]. In the former work, the authors stated that the average-intensity atlas was blurred around teeth owing to artifacts, which limited the capabilities of accurate non-rigid registration, and consequently segmentation, around dental structures. A single image Atlas, like proposed in this paper, will not suffer from these limitations. What we propose here is the use of the deformation field as the main feature to build an atlas. Atlas construction based on deformation field analysis has been an important topic of research [29]. More closely related to our work is the one in [11], where it is proposed to build an average shape and average intensity atlas from a set of subjects. In that work, one image is chosen as an initial reference from a set of five magnetic resonance images of the brain, and an iterative procedure based on deformation fields is applied to construct an image that would represent the shape and intensity average of the population. They have shown convergence for any choice of the reference image to the final average anatomic template. However, only five images were used for template construction; such limited database may not be enough to assess local shape differences or to generalize the fast convergence into an unbiased template. In this paper, contrary to [11], we do not perform intensity average. This is mainly motivated by the fact that the main modality in the context of head and neck radiotherapy is computed tomography (CT), in which dental implants may cause metal artifacts in the volumes, making impracticable to average intensities as performed in most proposed methods [28]. Also, instead of performing an iterative construction of the shape average, we propose a criterion and a procedure to find, among the original images in a database, the one that is closest to this shape average and use this image as the atlas. Our results, obtained on a database of thirty one CT images, show that the performance of this “closest to average” atlas in automatic segmentation is as good as the performance of the average itself. Moreover, this atlas does not suffer from blurring introduced by averaging and interpolation procedures. The paper is organized as follows. In the next section we briefly describe the non-rigid registration method, called Morphons, used in this work. In Section 3, we describe the atlas building methodology, as well the implementation details. The numerical results are presented in Section 4, along with an illustrative qualitative evaluation. An analysis of these results and their implications is given in Section 5. Finally, Section 6 presents concluding remarks, also pointing out future work. 2. Image registration
are CT volumes, the data tend to be quite noisy, containing artifacts that cause intensity variations. These artifacts are illustrated in Fig. 1, which shows slices of an axial and a sagittal view of one particular subject of our database. 2.1. The Morphons method The Morphons method is a non-rigid registration method using an iterative deformation scheme where the displacements estimations are found from local phase difference [30,14]. To find local phase information, a set of quadrature filters, each one sensitive to structures in a certain direction, is applied to both the fixed image and the atlas. The information on local phase allows the detection of local structures, such as lines and edges [7]. In order to detect large structures a multiresolution scheme is used. This multiresolution scheme involves an iterative accumulation of the deformation field under the influence of certainty measures, which are associated with the displacement estimates found in each iteration. The formulation of the Morphons method is detailed in the sequel. Let x = (x, y, z) denote a point in the space, I(x) denote an image and Q(x) denote the impulse response of a quadrature filter. When the quadrature filter is applied to an image, the output is given by q(x) = Q (x) × I(x) = A(x)ejϕ(x)
Since the impulse response of a quadrature filters is complex, the output has a magnitude, A(x), and a phase value, ϕ(x). Applying the quadrature filter to a pair of images IF (x) and IM (x) – representing the fixed and moving images respectively – results in q1 (x) = Q (x) × IF (x) = A1 (x)ejϕ1 (x) q2 (x) = Q (x) × IM (x) = A2 (x)ejϕ2 (x)
(2)
Now, compute the product of q1 (x) and the complex conjugate of q2 (x): q1 (x)q∗2 (x) = A1 (x)A2 (x)ejϕ(x)
(3)
where the phase difference is the argument of this product: ϕ(x) = ϕ1 (x) − ϕ2 (x). The Morphons method determines the local displacement estimate dı to be performed in a certain filter direction ı as being proportional to the local phase difference of the filter responses in that direction. That is, at each voxel dı ∝ ϕ(x). At each voxel, a displacement estimate is found for each filter in the filter set, that is, a displacement field dı is obtained for each ˆ ı . These fields are combined into one displacement filter direction n field by solving a least squares problem: min
d
Registration is the process of finding a transformation that best matches two images according to a measure of similarity. One image is the reference, which remains fixed during the registration process, while the other is deformed in the geometric space of the reference. The reference image is also called fixed or target image and the transformed image is called the moving image. For the atlas building, it is important to use very accurate registration methods. The non-rigid algorithm used in this application is the Morphons algorithm. This method was chosen based on our previous work [23] and also in [12], in which we have shown with quantitative results the capability of this non-rigid multi-modality registration method to segment regions of interest (ROI) in head and neck CT images. Morphons extracts special features from the intensity images (output of quadrature phase filters) and works on these features instead of the original intensity. This makes the Morphons algorithm appropriate for multi-modality registration and noisy applications. Even though in our application all images
(1)
T
ˆ ı d − dı )]2 [cı (n
(4)
ı
ˆ ı is the direction of filter where d is the sought displacement field, n ı, and cı is a certainty measure associated to the i-th direction. This certainty measure is usually taken as the magnitude of Eq. (3), that is, cı = A1 (x)A2 (x) [14]. Iterative accumulation of a dense deformation field is done under the influence of the certainty measures: da =
ca da + ck (da + dk ) ca + ck
(5)
where da indicates the updated accumulated deformation field, da is the accumulated field from the previous iteration and dk is the displacement estimate derived in the current iteration k. ca and ck are certainty estimates associated with the accumulated deformation field and the displacement estimates, respectively. ck is the sum of certainty values cı in Eq. (4) for all filter directions. ca is
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
449
Fig. 1. Slices of an axial view and a sagittal view of one image of the database showing metal artifacts from dental implants.
found by an accumulation procedure by using the certainty values as certainties of themselves: ca =
ca2 + ck2 ca + ck
(6)
This iteration is applied at each level of the multiresolution scheme. From one level to another, the deformation field is upsampled and accumulated until the final level is reached [30]. 3. Anatomical atlas The aim of atlas construction is to have a unique image/volume that represents a population’s anatomy. The main purpose of the atlas in head and neck radiotherapy is to automatically segment the zones at risk and also the regions with high probability of tumoral propagation in the fat tissues [10]. The union of all volumes of interests, including some error margin, will be designated in the following as clinical target volumes (CTVs). Several methods have been proposed to create a brain atlas, but the suitability of such atlases in other anatomical regions has neither been demonstrated nor validated. In fact, in our application to head and neck we encountered difficulties that do not exist in brain MR images; these difficulties severely limit the application of most proposed methods in head and neck CT. Firstly it is prohibitive to have CT images/volumes of normal subjects due to the nature of the exam. Another major problem is the presence of artifacts, which invalidates the use of intensity averaging techniques. Recently, an atlas methodology for constructing head and neck atlases was investigated using eight patients [3], based on [11], using three different non-rigid registration algorithms. The drawback of their work is the degree of blurriness presented in the final atlas [2]. 3.1. Atlas model Our atlas building consists in finding a subject that best represents anatomically a population. In order to chose this subject, one needs an unambiguous numerical criterion which indicates that a certain image is the most representative in size/shape of human anatomy, this image being then chosen as the atlas. Non-rigid registration accounts for the differences of subjects due to size/shape variation, enabling the deformation of each of the subjects to be quantitatively and qualitatively compared. Accordingly, we will use deformation fields as the main feature to discriminate the best subject to be the atlas. The subject for which the average of the deformation field is closest to zero, after having all other subjects being registered into him, is the optimal shape
template or atlas. This is the principle of our method, assuming that all the images in the database are geometrically aligned. Let us now introduce formally our methodology. A deformation field is a function that represents a correspondence between two different subjects. Let D(x) = (Dx (x), Dy (x), Dz (x)) represent a deformation field for the whole image/volume, Dx (x), Dy (x) and Dz (x) being its components in each spatial direction. Let i and j be the indices of two subjects in a database of N images. The displacements from subject j to subject i in each direction x, y and z are described respectively by Dxij (x), Dyij (x) and Dzij (x). This triplet defines the deformation field that registers subject j into subject i and represents how much a subject has to be deformed to have the same size/shape of the other. After all subjects j have been registered in a certain subject i, the average of all the deformation fields for subject i will be found. The average displacements in each direction x, y and z are, respectively:
Dxi (x) =
Dyi (x) =
Dzi (x) =
1 N−1
1 N−1
1 N−1
N
Dxij (x)
(7)
Dyij (x)
(8)
Dzij (x)
(9)
j=1,j = / i N j=1,j = / i N j=1,j = / i
and Di (x) = (Dxi (x), Dyi (x), Dzi (x)) is the average deformation field for the subject i. The subject i whose average of the deformation field is closest to zero is the optimal shape template. To implement this rationale, we need to measure the closeness to zero in a concise manner, since Di (x) is a very dense information. We do that by using the norm of the magnitude, as defined in the following equations.
|Di (x)| =
Dxi (x)2 + Dyi (x)2 + Dzi (x)2
1 |Di (x)| m×n×p m
Di =
n
(10)
p
(11)
x=1 y=1 z=1
where m, n and p are the image dimensions. The quantity Di thus defined represents the mean bending energy of the deformation field for subject i. It represents, in a single number, how far the subject i stands from the average anatomy. Hence,
450
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455 Table 1 D values (in voxels) for each patient of the database.
Fig. 2. Scheme that illustrates the atlas building process. See text for details.
i
Subject
Di
i
Subject
Di
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 31
Subject1 Subject2 Subject3 Subject4 Subject5 Subject6 Subject7 Subject8 Subject9 Subject10 Subject11 Subject12 Subject13 Subject14 Subject15 Subject31
1.79 3.59 1.07 2.12 2.62 2.09 2.28 4.45 2.23 1.41 2.12 2.16 1.42 2.94 1.32 5.31
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Subject16 Subject17 Subject18 Subject19 Subject20 Subject21 Subject22 Subject23 Subject24 Subject25 Subject26 Subject27 Subject28 Subject29 Subject30
4.32 2.36 2.51 1.61 1.96 1.72 1.48 4.28 1.87 1.80 1.65 2.24 1.62 1.49 1.50
the subject that requires the minimum bending energy is the optimal anatomical atlas A, that is: A = Ik ,
where
k = arg(minDi )
(12)
4.1. Database
i
The block diagram in Fig. 2 illustrates the atlas building process. In this scheme, let each gray box be a subject in a database of 6 images. The boxes are disposed by their size/shape in relation to the central image, this latter being the most central anatomy in this database. Since we do not know a priori which image has the most representative size/shape, each and everyone at a time will be the fixed image and all the 5 subjects left will be non-rigidly registered into him. Based on the left diagram of Fig. 2, subject S1 starts being the atlas candidate. The deformation field resulting from registering image S2 into S1 is denoted by D12 , registering image S3 into S1 yields the deformation field D13 , and so on. After all images have been registered to S1, Eqs. (7), (8) and (9) are calculated for the set of deformation fields D1j from all registrations pairs 1j. Then the Eq. (11) is computed for subject S1. Next, all steps described above have to be repeated considering S2 as the fixed image – or the atlas candidate – as shown at the right in Fig. 2. The procedure is repeated, each subject at a time being the atlas candidate. The subject i that yields the minimum value, as specified in Eq. (12), is the best shape template for that population and is the proposed atlas A. In order to illustrate our atlas construction procedure, we took a real image of a hand and artificially created two others to have an enlarged and a reduced version of it. Fig. 3 shows the real hand on top, the enlarged version in the left bottom and the reduced one in the right bottom. In this experience, after registering both modified hands at the original one, we calculated the mean bending energy of the deformation field, finding a value of D1 = 0.44. Registering the enlarged hand and the original one into the reduced hand, the mean bending energy provided a value of D2 = 8.42. Then the reduced hand and the original one were registered into the enlarged hand, providing a mean bending energy value of D3 = 6.48. The mean bending energy of the average deformation field was around 15–20 times smaller considering the original hand as the central anatomy. This corroborates the expected behavior of our method, showing that our metric is capable of finding the image with a central anatomy in the population.
4. Results We describe in the sequel the application of the proposed atlas building procedure to a real database of head and neck images, as well as the atlas’ use in automatic segmentation of OARs and CTVs.
A database of 31 patients of 3D CT images previously segmented by a radiologist has been used to build the atlas. The size of the CT volumes is of 256 × 256 × 128 voxels with a voxel size of 0.9765 mm × 0.9765 mm × 2.1093 mm. All images in this database are male patients with age between 50 and 70 years old.
4.2. Implementation Before performing Morphons registration, all images in the database have to be at same geometric space; so we firstly apply a rigid registration. The geometric alignment is obtained by minimizing the mutual information between the fixed and the moving image [20] using a simultaneous perturbation stochastic approximation (SPSA) [4]. The rigid transformation is a superposition of a 3-D rotation and a 3-D translation and the registration parameter is a six-component vector consisting of three rotation angles and three translation distances [18]. The rigid registration is done only once, putting all database at the same geometric space of any chosen subject. The choice of this subject can be arbitrary, because linear transformation such as the rigid registration preserves the anatomy of the moving image and consequently does not affect the result. The rigid registration algorithm was performed using the ITK environment (http://www.itk.org) and Morphons was performed using Matlab. A total of 31 rigid registrations and 31×30 non-rigid registrations has been performed. For each image i being the fixed in the registration process, we have calculated D from Eq. (11). Table 1 summarizes the numerical results – the D values for all 31 patients, which are given in unities of voxels. It is seen that subject 3 was found to give the minimal value of D = 1.07. This means that subject 3 requires the minimum mean bending energy for all database to be deformed into him, and that the average displacement for each voxel is a distance equivalent to 1.07 voxels. Fig. 4 shows an axial and sagittal view of subject 3. To illustrate that subject 3 has indeed a central anatomy regarding the size/shape in relation to others in the database, we took two subjects with high values of D: patient 31 with D31 = 5.31 and patient 16 with D16 = 4.32. These high values mean that these subjects have extreme anatomies, that is, anatomies that are significantly distant from the population’s average. This is indeed the case, as we can observe in the Fig. 5, which presents an axial and a sagittal view of patient 31, and in Fig. 6 which presents an axial and a sagittal view of patient 16. This means that in average, those subjects require higher mean bending energy.
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
Fig. 3. Database of hands. Top: original hand; left bottom: artificially enlarged hand; right bottom: artificially reduced hand.
Fig. 4. Slice of patient 3. Left: axial view; right: sagittal view.
Fig. 5. Slice of patient 31. Left: axial view; right: sagittal view.
451
452
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
Fig. 6. Slice of patient 16. Left: axial view; right: sagittal view.
4.3. Atlas based segmentation Next, we apply our atlas to the task of automatic segmentation of the regions of interest in head and neck radiotherapy and compare the results to the manual segmentation. The regions of interest were manually contoured in all images by a physician, including the right and left parotid glands, the spinal cord, the clinical target volume and also the body contour. From each contour a binary image (mask) has been created, setting ones inside the contour (ROI) and zeroing the remainder of the grayscale image. We use this manual segmentation performed by an expert for comparison purposes. The atlas based segmentation requires that the atlas itself has been previously segmented, so that binary images are available for each region of interest. The segmentation of a given subject’s image is then achieved in two steps. First the atlas is registered into the subject and then the resulting deformation field is applied to the binary masks of the atlas. The result of this process is a set of binary masks for the patient, each one representing one region of interest. To evaluate the segmentation done by the atlas, the deformed masks are compared to the manual segmentation done by the physician. The metric used to measure spatial intersection of two masks was the similarity index (SI) [31]: SI = 2 ·
|MP ∩ MA | |MP | + |MA |
(13)
where | · | is the number of non-zero voxels of the masks and SI ∈ [0, 1]. The SI was calculated between the masks originated from the manual segmentation – MP – and the masks originated from the automatic segmentation using our atlas – MA – as in Eq. (13). Table 2 presents the resulting SIs for the ROIs of selected subjects from the database.1 The procedure in [11] constructs iteratively an average shape of the population, whereas we only pick, among the available (real) images, the one that is closest to this average. To compare the atlas obtained with our strategy with the one described in their work, we unbiased our atlas using their method. This is achieved by applying the inverse average deformation field into the first image chosen as the atlas, creating a new atlas. Then all images in the database are registered into this new version of the atlas, and then again the inverse average deformation field is applied to the present atlas, creating again a new version of the atlas. This process is done iteratively until convergence is reached, that is, until the new atlas does not change significantly from the previous iteration. This is also similar to the method proposed in [15]. This scheme is illustrated in Fig. 7. In this figure, the symbol A indicates the atlas A proposed here and the boxes represent the
1 The non-filled lacunas in Table 2 correspond to the contours that were not available in the Dicom-RT original images.
other images in the database. At the first iteration, all subjects, denoted by S1 . . . S6 in Fig. 7, are registered into the atlas A. Then, the average deformation field is calculated and its inverse is applied to the atlas A, creating a new version of an atlas, denoted by G1 . These steps are repeated, but now registering all images from the database in G1 , creating a new estimation of the average (Guimond’s) atlas, G2 . The convergence of our atlas A to the Guimond’s atlas occurred after 6 iterations and we called it atlas G. To illustrate it qualitatively, an axial and a sagittal views are shown of Guimond’s atlas in Fig. 8. To help the visualization of the atlas G, the A is also presented in Fig. 9 with a blue contour placed at the bones, and this contour is overlayed in G to facilitate the visualization of anatomic differences between those images. After atlas G was created, it was used to segment the ROIs previously described. The results are summarized in Table 3, where the same SI metric was used. It can be seen in this table that the performance of both atlases in the automatic segmentation task is quite similar. Indeed, performing an analysis of variance (ANOVA) to these data with a level of significance of 5% established that the results obtained with the two different atlases present no statistical differences (p-values p 0.05 were found). 5. Discussion In reference [11] an average atlas for the brain is arrived at after three iterations. However they have used a small database of only five patients for building the average shape. A similar work has reported more iterations for a larger database [13]. In this paper, we achieved convergence to the average atlas G after six iterations, with the major differences that our database consists of thirty one subjects and the first image used to build the Guimond’s atlas was not a random image from our database. The size of the database regarding the population that the atlas is supposed to represent is an important issue. In this work we have used a reasonably large database (thirty one subjects) to represent a somewhat restricted population, namely only male adults, so we expect the database to be sufficiently representative. The mean bending energy of our atlas A is 1.07, meaning that each voxel must be displaced an average of 1.07 voxels when registering all the subjects of the database in the atlas. Applying an iterative procedure to search for the average anatomy results in an alternative atlas whose final mean bending energy, after six iterations, is 0.20. Although this resulting atlas is certainly close to the central anatomy of the population, this difference does not seem to be relevant. Indeed, it has been seen that the performance of the two atlases in automatic segmentation is not distinguishable, and that anatomical differences between them are hardly perceivable for the naked eye. Hence our results imply that, in the quest for the average anatomy, there is a limiar beyond which it is useless
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
453
Table 2 Similarity index for segmenting the regions of interest from head and neck images using the proposed atlas A versus manual delineation done by a physician. ROI
Pat4
Pat10
Pat13
Pat16
Pat17
Pat23
Pat29
Body Spinal cord Right parotid Left parotid Right CTV Left CTV
0.98 0.69 0.70 0.68 – 0.73
0.96 0.74 0.69 0.57 – 0.47
0.97 0.65 0.54 0.30 0.60 0.57
0.97 0.55 0.70 – 0.70 0.64
0.96 0.51 0.54 0.49 0.55 0.59
0.98 0.77 0.63 0.65 – 0.70
0.98 0.73 0.60 0.57 0.65 –
Fig. 7. Scheme illustrating 3 iterations of Guimond’s method in the proposed Atlas A, generating an unbiased Atlas G.
Fig. 8. Axial and sagittal views of the atlas Guimond at the sixth (final) iteration [11]. The same blue contour from Fig. 9 was overlayed this image to help the visualization of anatomic difference between the two images. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Fig. 9. Axial and sagittal views of the atlas A, also used as the reference image for constructing Guimond [11]. A blue contour was placed in the bone structure to help visualization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Table 3 Similarity index for segmenting the regions of interest of head and neck images using atlas G versus manual delineation done by a physician. ROI
Pat4
Pat10
Pat13
Pat16
Pat17
Pat23
Pat29
Body Spinal cord Right parotid Left parotid Right CTV Left CTV
0.97 0.68 0.66 0.60 – 0.70
0.96 0.74 0.68 0.54 – 0.47
0.97 0.64 0.54 0.31 0.56 0.56
0.96 0.54 0.68 – 0.64 0.54
0.96 0.52 0.56 0.48 0.57 0.57
0.97 0.73 0.55 0.63 – 0.66
0.97 0.70 0.60 0.60 0.68 –
454
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455
to proceed. On the other hand, most methods of atlas building use some sort of averaging and/or interpolation, which induces blurring; it is thus advisable not to proceed any further than necessary in such methods. In our database this limiar is slightly above 1 voxel and any image whose mean bending energy is smaller than this is very likely to be a good atlas. Whether this particular value can be taken as a general rule is a topic of further investigation. Recent results strongly indicate that the strategy of multi-atlas based segmentation is more effective than a single atlas to cope with the anatomy variability [9,24]. In the former work, clusters are created and then an average image of each cluster is an atlas. The methodology for atlas building proposed here could be used to create the cluster based on the mean bending energy for a multiatlas strategy. Instead of creating an average image, each cluster would have a median image chosen by our criterion that would be neither biased nor blurred. Also, a new image of a patient to be segmented would select the closest atlas based on the least mean bending energy. Our atlas could also be used as the initialization for a more sophisticated method of segmentation such as the active contour based segmentation as in [8], in which they have chosen an arbitrary image to be the atlas. 6. Conclusions and final remarks In this paper we have proposed a methodology to build an atlas that represents the size and shape variability present in human anatomy based on dense deformation fields using Morphons. In our previous work we have shown that Morphons is an effective strategy for performing non-rigid registration of the head and neck anatomy. Patients with high mean bending energy – a large D value – have more difference in anatomy in relation to the others, as we have shown. The atlas proposed here consists of the subject with the lowest D value within a population. The main advantage of our methodology is that the chosen atlas is an original image/volume and thus does not have any degree of blurriness imposed by interpolation or iterative processes, which could impose problems in future registration when using it for atlas based segmentation. On the other hand, it may happen, particularly when a restricted database is considered, that even the lowest D is still not small enough – meaning that the subject closest to the average anatomy is not close enough to it in order to be a good atlas. In such cases our atlas could be used as a complement to iterative procedures, by providing the best initial condition for them. One of our future work objectives is to improve the SI for some regions of interest being investigated. The goal is to build a statistical atlas based on [27] and also investigate the use of segmentation methods based on active contours [8], which the atlas based segmentations will be used as a initialization. Alternative metrics, such as in [1], will also be considered in future works. Finally we need to mention the computational complexity of our method. The number of registrations to be computed is N · (N − 1). Although this method is computationally intensive this does not pose any serious problems, since our calculations are performed off-line and all the registrations has to be done only once. Acknowledgements We gratefully acknowledge the support from the SIMILAR Network of Excellence, the Brazilian Education Ministry (MEC-CAPES) and the National Council of Technological and Scientific Development (CNPq) for financial support. The authors are grateful to Dr. Vincent Grégoire, Xavier Geets and Pierre Castadot for providing the images and their manual delineations. The authors wish to thank Vincent Nicolas for developing efficient visualization tools (http://www.medicalstudio.org/).
References [1] V. Arsigny, O. Commowick, X. Pennec, N. Ayache, A log-Euclidean framework for statistics on diffeomorphisms, in: Proceedings of the 9th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI’06), 2–4 October, Part I. No. 4190, 2006. [2] O. Commowick, V. Grégoire, G. Malandain, Atlas-based delineation of lymph node levels in head and neck computed tomography images, Radiother. Oncol. 87 (2) (2008) 281–289. [3] O. Commowick, G. Malandain, Evaluation of atlas construction strategies in the context of radiotherapy planning., in: Proceedings of the SA2PM Workshop (From Statistical Atlases to Personalized Models), Copenhagen, held in conjunction with MICCAI 2006, October 2006. [4] M.D. Craene, Dense deformation field estimation for pairwise and multisubjects registration, Ph.D. thesis, Université Catholique de Louvain, B - 1348 Louvain-la-Neuve Belgique, August 2005. [5] M.B. Cuadra, C. Pollo, A. Bardera, O. Cuisenaire, J.-G. Villemure, J.-P. Thiran, Atlas-based segmentation of pathological MR brain images using a model of lesion growth, IEEE Trans. Med. Imaging 23 (10) (2004) 1301–1314. [6] B.M. Dawant, S.L. Hartmann, J.-P. Thirion, F. Maes, D. Vandermeulen, P. Demaerel, Automatic 3D segmentation of internal structures on the head in MR images using a combination of similarity and free form transormations: Part I, metholody and validation on normal subjects, IEEE Trans. Med. Imaging 18 (10) (1999) 909–916. [7] W.T. Freeman, E.H. Adelson, The design and use of steerable filters, IEEE Trans. Pattern Anal. Mach. Intell. 13 (9) (1991) 891–906. [8] S. Gorthi, V. Duay, N. Houhou, M. Bach Cuadra, U. Schick, M. Becker, A.S. Allal, J.P. Thiran, Segmentation of head and neck lymph node regions for radiotherapy planning, using active contour based atlas registration, IEEE J. Sel. Top. Signal Process. 3 (2009) 135–147. [9] S. Gorthi, M. Bach Cuadra, U. Schick, P.-A. Tercier, A.S. Allal, J.-P. Thiran, Active multi-atlas based segmentation of head and neck CT images using active contour framework, in: MICCAI workshop on 3D Segmentation Challenge for Clinical Applications, 2011, p. 2010. [10] V. Grégoire, P. Levendag, J. Bernier, M. Braaksma, V. Budach, C. Chao, E. Coche, J. Cooper, G. Cosnard, A. Eisbruch, S. El-sayed, B. Emami, C. Grau, M. Hamoir, N. Lee, P. Maingon, K. Muller, H. Reychler, CT-based delineation of lymph node levels and related CTVs in the node-negative neck: DAHANCA, EORTC, GORTEC, NCIC,RTOG consensus guidelines, Radiother. Oncol. 69 (3) (2003) 227–236. [11] A. Guimond, J. Meunier, J.-P. Thirion, Average brain models: a convergence study, Comput. Vision Image Understanding 77 (2) (2000) 192–210. [12] G. Janssens, J. Orban de Xivry, B. Macq, Evaluation of nonrigid registration models for interfraction dose accumulation in radiotherapy, Med. Phys. 36 (9) (2009) 4268–4276. [13] S. Joshi, B. Davis, M. Jomier, G. Gerig, Unbiased diffeomorphic atlas construction for computational anatomy, Neuroimage 23 (Suppl. 1) (2009). [14] H. Knutsson, M. Andersson, Morphons: segmentation using elastic canvas and paint on priors, in: IEEE International Conference on Image Processing (Proceedings of the IEEE International Conference on Image Processing’05), Genova, Italy, September 2005. [15] P. Kochunov, J.L. Lancaster, P. Thompson, R. Woods, J. Mazziotta, J. Hardies, P. Fox, Regional spatial normalization: toward an optimal target, J. Comput. Assist. Tomogr. 25 (5) (2001) 805–816. [16] P. Lorenzen, M. Prastawa, B. Davis, G. Gerig, E. Bullitt, S. Joshi, Multi-modal image set registration and atlas formation, Med. Image Anal. 10 (June (3)) (2006) 440–451. [17] M.D. Craene, M.D. Duay, V. Macq, B. Pollo, C.J. Thiran, Dense deformation field estimation for atlas-based segmentation of pathological MR brain images, Comput. Methods Programs Biomed. (August) (2006) 66–75. [18] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Trans. Med. Imaging 16 (2) (1997) 187–198. [19] P. Maingon, L. Ramus, G. Odin, V. Grégoire, V. Darcourt, P.Y. Marcy, N. Guevara, M.H. Orlanducci, S. Marcie, G. Poissonnet, A. Bozec, O. Dassonville, F. Demard, L. Castillo, J. Santini, G. Malandain, DENTALMAPS: Automatic Dental Delineation for Radiotherapy Planning in Head and Neck Cancer, Int. J. Radiat. Oncol. Biol. Phys. (2011). [20] D. Mattes, D.R. Haynor, H. Vesselle, T.K. Lewellen, W. Eubank, Pet-ct image registration in the chest using free-form deformations, IEEE Trans. Med. Imaging 22 (1) (2003) 120–128. [22] H. Park, P. Bland, C.R. Meyer, Construction of an abdominal probabilistic atlas and its application in segmentation, IEEE Trans. Med. Imaging 22 (4) (2003) 483–492. [23] A. Parraga, A. Susin, J. Pettersson, B. Macq, M.D. Craene, Quality assessment of non-rigid registration methods for atlas-based segmentation in head-neck radiotherapy, in: Proceedings of IEEE International Conference on Acoustics, Speech, & Signal Processing, Honolulu, Hawaii, USA, April 2007. [24] L. Ramus, Conception et utilisation d’atlas anatomiques pour la segmentation automatique: application à la radiothérapie des cancers ORL, Ph.D. thesis, Université de Nice Sophia Antipolis, July 2011. [25] A. Rao, R. Chandrashekara, G.I. Sanchez-Ortiz, R. Mohiaddin, P. Aljabar, J.V. Hajnal, B.K. Puri, D. Rueckert, Spatial transformation of motion and deformation fields using nonrigid registration, IEEE Trans. Med. Imaging 23 (9) (2004) 1065–1076.
A. Parraga et al. / Biomedical Signal Processing and Control 7 (2012) 447–455 [27] D. Rueckert, A.F. Frangi, J.A. Schnabel, Automatic construction of 3D statistical deformation models of the brain using non-rigid registration, IEEE Trans. Med. Imaging 22 (8) (2003) 1014–1025. [28] D. Seghers, E. D’Agostino, F. Maes, D. Vandermeulen, P. Suetens, Construction of a brain template from mr images using state-of-the-art registration and segmentation techniques., in: MICCAI (1), 2004, pp. 696–703. [29] Q. Wang, D. Seghers, E. D’Agostino, F. Maes, D. Vandermeulen, P. Suetens, A. Hammers, Construction and validation of mean shape atlas templates for
455
atlas-based brain image segmentation., in: Proceedings Information Processing in Medical Imaging – IPMI, vol. 3565, 2005, pp. 689–700. [30] A. Wrangsjö, J. Pettersson, H. Knutsson, Non-rigid registration using morphons., in: Proceedings of the 14th Scandinavian Conference on Image Analysis (SCIA’05), Joensuu, 2005 June. [31] A. Zijdenbos, B. Dawant, R. Margolin, A. Palmer, Morphometric analysis of white matter lesions in MR images: method and validation, IEEE Trans. Med. Imaging 13 (December (4)) (1994) 716–724.