In this work, we present multi-atlas based techniques for both segmentation and landmark detection. We focus on modality and anatomy independent techniques to be applied to a wide range of input images, in contrast to methods customized to a specific anatomy or image modality. For segmentation, we use label propagation from several atlases to a target image via a Markov random field (MRF) based registration method, followed by label fusion by majority voting weighted by local cross-correlations. For landmark localization, we use a consensus based fusion of location estimates from several atlases identified by a template-matching approach. Results in IEEE ISBI 2014 VISCERAL challenge as well as VISCERAL Anatomy1 challenge are presented herein.
For segmentation, a multi-atlas based technique is used by registering several atlases individually to a target image using our implementation of the MRF-based deformable registration method in [GKT+08]. These registrations are then used to propagate the anatomical labels (ground-truth annotations) from each atlas image into the target coordinate frame. At voxel level, a majority voting is held to decide the winning label where each label votes based on the locally-normalized cross-correlation (LNCC) [CBD+03] of the registered atlas to the given target at that location. 2.1
Finding an optimal displacement vector field Tˆ can be defined as the minimization of a functional: Tˆ = arg min E (T, X, An)
T SX,n = SAn (Tˆ).
Although MRF-based registration is a relatively robust method, it can only guarantee a locally optimal solution and is therefore susceptible to poor initialization. Furthermore, for a grossly different atlas, correspondences for registration may not be guaranteed. Accordingly, the segmentation from where X is a target image and An is an atlas image and E is the registration energy. MRFs provide an efficient means for the minimization of such energy, with the main advantage being that it does not rely on the gradient of the criterion and therefore is less prone to poor local optima. In order to use MRFs for solving the minimization problem, this energy is decomposed into unary (ψ) and pairwise (Ψ) potentials over discrete labels as follows: (1) (2) (3) E (T, X, An) = X p∈Ω ψp (lp) +
X q∈N (p)
λ Ψpq (lp, lq) , where Ω is the discretized image space. The continuous displacement space is sampled discretely, so that each registration label lp and lq in the set of all registration labels LR maps to a unique displacement vector d~p. The unary potentials are then a local similarity metric, measuring the fit between the deformed atlas image and the target image at a location i. The pairwise potentials correspond to prior assumptions over the displacement, which is often implemented as a smoothing over the neighborhood N as justified by a first-order Markov assumption. λ is the pairwise weight.
For robust and smooth solution of (2), an efficient method was proposed in [GKT+08] that seeks the displacements d~ of control points in a multi-resolution cubic B-spline framework. We use our implementation of this method with four levels of detail, where the coarsest grid resolution has three nodes along the shortest edge of an input image and a spacing as isotropic as possible given that a control-point is required on each corner of the image. Each following level of the resolution hierarchy has twice the resolution compared to the previous step. At each level of detail, we sample four displacements in each cardinal direction, yielding 25 displacement samples in total. Within each direction, samples are equidistant, with the largest displacement set to 0.4 times the control grid spacing. This was shown to guarantee diffeomorphic deformations [RAH+06]. For each level of detail, we re-run the MRF registration with the displacements being re-scaled by the golden ratio 0.618. The resulting displacements are then composed onto the previous deformation. This guarantees that the result is still diffeomorphic and sub-pixel accuracy can be achieved. For the unary potentials we use normalized cross correlation (NCC) of patches centered around each control point with their radius equivalent to control grid spacing. Euclidean distance between displacements of neighboring control grids is used to penalize non-smooth deformations. Tree-reweighted message passing (TRW-S) [Kol06] is employed to find a solution to each energy minimization instance.
A segmentation candidate of the target image X based on the atlas An is then obtained by applying the resulting displacement field Tˆ to the known segmentation SAn as follows: a single atlas may not be satisfactory. It was shown in different research fields that combination of multiple weak information sources can surpass average accuracy. Multiple segmentations from different atlases were combined in [HHA+06].
Assume that N segmentation candidates of a target image X are computed from N atlases via (3). Let final target segmentation SX be an image of the same size as X, where pixels take on values from the set LS={1, . . . , NS} such that each discrete value corresponds to an organ or anatomical structure. An intuitive and straight-forward method to combine multiple segmentation estimates is then to choose the most frequent segmentation label (majority voting, MV) at each location p: SXMV(p) = arg max X δ (lS, SX,n(p)) .
lS∈LS n Such majority voting does not take into account the individual quality of each registration and therefore the resulting segmentation. We assume that post-registration image similarity between the deformed atlas and the target image is an indicator of segmentation reliabilty and can be used to locally assign weights w to each individual segmentation. The resulting weighted majority vote (wMV) can then be formalized as follows:
SXwMV(p) = arg max X wn(p) δ (lS, SX,n(p)) .
lS∈LS n To obtain the weights w, we use local normalized cross correlation (LNCC,[CBD+03]) between image X and deformed atlas An(Tn). The advantages of LNCC are its smoothness and fast computation time due to convolution with Gaussian kernels:
LNCC(X, Y, p) =
hX, Y i(p) σX (p) σY (p) X = GσG ∗ X hX, Y i(p) = X · Y (p) − X(p) · Y (p) σX2 (p) = X2(p) − X2(p), where ∗ is the convolution operator and GσG is a Gaussian kernel with standard deviation σG . From the LNCC metric, we compute the weights: wn(p) = 1 − LNCCσ(X, An(Tn), p) γ 2 (4) (5) (6) (7) which normalizes LNCC to the range [0, 1]. γ is used to scale the similarity such that contributions from individual segmentations are well spread [IK09]. 3
For anatomical landmark detection, we use a template based approach from multiple atlases, the location estimates from which are fused based on their consensus. We localize each landmark ℓ separately from the others using the two stages below. To localize the unknown voxel coordinates pℓ of landmark ℓ in the target image X, we perform the following template matching procedure from each atlas An where n represents the atlas index. 3.1
The template is set as a box-shaped image region Aℓn in the current atlas. Similarly, a box-shaped search region Xℓ is defined in the target image. Both such regions are chosen targeting a physically isotropic region of interest (ROI) in corresponding image, while limiting the maximum number of ROI voxels to ensure efficient computation. Specifics of ROI selection are given in Table 1. pℓAn · ||AXn||
When setting the template size, the trade-off between it containing sufficient image features and final localization precision was considered. Template half-width was set to 20 mm empirically via cross-validation in multiple modalities using different template sizes (e.g., 10, 20, 30, ...). The search region is centered around a gross estimate of the landmark location, which is the normalized voxel coordinates of the landmark from the atlas; using the fact that both the atlas and the target have similar fields of view (ie. both abdomen, thorax, or whole body). Note that our large search region covers most or all of the image in many modalities (e.g. in MRce) or at least a quadrant thereof (e.g. in CT), such that the searched landmark can be guaranteed to exist therein. For template matching, the template is convolved over the search region by computing two independent similarity metrics, sum of squared differences (SSD) and normalized cross-correlation (NCC), at each template location i with respect to search image. Both values are then normalized linearly to [0, 1] such that they are both 1 at the best match location. A combined similarity metric SSDa·CORb is then computed, where the parameters a=2 and b=3 were determined empirically via cross-validation with several powers. The maximum of this combined metric gives the best match location estimate pℓn for landmark ℓ considered atlas An. 3.3
From cross-validation trials with different techniques such as the mean and weighted average of location estimates, the median operator was determined to be the best method for fusing location estimates. Accordingly, each axis coordinate of the the target landmark location is found as the median value of those axes from atlas estimates. The entire process can be summarized as: • Repeat for each atlas An : → Crop landmark template Aℓn centered at given landmark location pℓAn in atlas image An → Crop a large search region Xℓ centered around a grossly approximated location in X → Compute SSD Aℓn, Xℓ and COR Aℓn, Xℓ → pℓn = arg max SSD[i]2 · COR[i]3
i • pℓ = median pℓn | ∀n 4
Throughout the results, the following abbreviations are used for the image modalities: whole-body CT images (CT), thorax+abdomen contrast-agent CT images (CTce), abdominal T1-weighted contrast-agent MR images (MRce), and whole-body T1-weighted MR images (MR).
We treated each modality separately. We used the training images from Anatomy1 benchmark as atlases, ie. depending on the modality of the test image, six atlases for CTce and seven for CT, MR, MRce. For convenience, we combined all organ segmentations for each atlas into a single multi-label segmentation image, which was then deformed using the atlas-to-target registration Tˆ described in Sec. 2.1. In Tab. 2, Dice overlap metric results regarding our segmentation approach reported by VISCERAL for the test images can be seen both for the VISCERAL Anatomy1 benchmark and the ISBI challenge.
In order to compare our technique to other participants’ in the ISBI challenge, we computed an average rank per organ per participant. For given test image, we assigned a rank to each method, ie. {1,2,...,P } where P is the number of participants submitted a (non-blank) output for that test image. In Fig. 1 the average of such ranks for all given test images is seen per anatomy. Our submission was the only entrant that aimed to segment all images in all modalities and it achieved competitive results for many organs as seen in the given figure. We did not plot ranks for the MR modalities, since we were the only participant to submit such results.
It is notable that the multi-atlas fusion has significantly lower segmentation accuracy for organs which has low volume, e.g. urinary, gall bladder, and the adrenal glands. One possible explanation is that such small structures are difficult to register using full-body images. Due to such misalignments, overlap between multiple deformed atlas segmentations for such label can be small, resulting in the weighted majority voting not selecting that label. participant2 participant3 participant4 our Approach participant1 participant2 participant3 participant4 our Approach does not require any customisation to a specific modality or organ competed in all categories and ranked satisfactorily compared to the results of other participants. [CBD+03] Pascal Cachier, Eric Bardinet, Didier Dormont, Xavier Pennec, and Nicholas Ayache.
Iconic feature based nonrigid registration: the PASHA algorithm. Computer Vision and Image Understanding, 89(2-3):272–298, February 2003. [GKT+08] Ben Glocker, Nikos Komodakis, Georgios Tziritas, Nassir Navab, and Nikos Paragios.
Image Analysis, 12(6):731–741, December 2008. Medical Dense image registration through MRFs and efficient linear programming. [HHA+06] Rolf a Heckemann, Joseph V Hajnal, Paul Aljabar, Daniel Rueckert, and Alexander Hammers. Automatic anatomical brain MRI segmentation combining label propagation and decision fusion. NeuroImage, 33(1):115–26, October 2006. [IK09] [Kol06]