Using the Logarithm of Odds to Define a Vector Space on Probabilistic Atlases

2.1 An Introduction to LogOdds

LogOdds are an example of a class of functions that map the space of discrete distributions (Kendall and Buckland, 1976) to Euclidean space. Let ℙ_M be the open probability simplex (the space of discrete distributions) for M labels:

Note that ℙ_M is an M−1 dimensional space as the M^th entry can be computed from the first M−1 entries. Furthermore, the space is open avoiding the degenerative distributions that are certain about the assignments. For the specific case of M = 2, ℙ₂ is the Bernoulli distribution ℙ ≜ {p|p ∈ (0, 1)} (Evans et al., 2000). Many binary classification problems use the Bernoulli distribution where p represents the probability that a voxel belongs to a particular anatomical structure and its complement p̄ = 1−p the probability of the voxel being in the background.

The multinomial LogOdds function logit(·) : ℙ_M → ℝ^M−1 of a discrete distribution p ∈ ℙ_M is defined as the logarithm of the ratio between the i-th and last entry of p:

with i ∈ {1,…, M − 1}. For the Bernoulli distribution, this function simplifies to the logarithm of the ratio between the probability p and its complement:

The inverse of the log odds function logit(·) is the generalized logistic function¹

where Z ≜ 1 + Σ_{j=1,…,M−1} ^et_j is the normalization factor.

Let 𝕃_M−1 be the M−1 dimensional space of LogOdds induced from ℙ_M:

We note that 𝕃_M−1 is equivalent to a (M−1)D real vector space and is thus a vector space.

2.2 The Relationship between LogOdds and PA

The function logit(·) and its inverse comprise a homeomorphism between ℙ_M and 𝕃_M−1 so that we can borrow the vector space structure on 𝕃_M−1 to induce one on ℙ_M.

3.1 Signed Distance Maps, LogOdds and Probabilistic Atlases

LogOdds maps in ${\mathbb{L}}_1^n$ define the boundary of a shape as a zero-level set function. One subset of maps in ${\mathbb{L}}_1^n$ are SDMs that also conform to the Eikonal equation (Rauch, 1991) with uniform speed. Thus, SDMs can always be interpreted as LogOdds maps, but the reverse is in general not true. The corresponding signed distance map transformation $\mathcal{D}(\cdot ):{\mathbb{B}}_2^n\to {\mathbb{L}}_1^n$ can be seen as a direct mapping between binary maps and the LogOdds space. Transformed to ${\mathbb{P}}_2^n$, these maps define probabilistic atlases, where voxels inside the object are represented by probabilities higher than 0.5 for the foreground and voxels outside with probabilities higher than 0.5 for the background. In the case of discrete data, the mapping ${𝒟}_{M-1}:{\mathbb{B}}_M^n\to {\mathbb{L}}_{M-1}^n$ is defined by combining a set of SDMs 𝒟 into a vector 𝒟_M−1(B) = (𝒟(B1), … ,𝒟(B_M−1)), where $B_j\in {\mathbb{B}}_2^n$ is the binary map corresponding to label j in B. This mapping is illustrated for a few label maps in the second row of Figure 3. Within the LogOdds space, SDMs not only represent shapes but also define a set of PAs in ${\mathbb{P}}_M^n$ (see middle row² of Figure 4), which are generated via the generalized logistic function 𝒫(·). These PAs characterize the inside of an object with higher probabilities than voxels outside the object.

The proper choice of PAs depends on the application. For example, one can generate PAs directly via Gaussian smoothing of each binary map B_j. This results in maps with values between zero and one. We then have to normalize the resulting maps to create a discrete distribution map, called GAUSS, that sums to one at each voxel location. An example of such a PA is shown in the last row of Figure 4, where we used a Gaussian filter of standard deviation 10. The corresponding LogOdds (last row of Figure 3) are generated via the logarithm of odds function logit(·).

In the case of Bernoulli distributions, these representations are very similar (first two columns of Figure 3 and Figure 4). However, once we turn to discrete distributions, stronger differences appear. The PAs generated from SDMs show distortions at the interface between the two circles (see second row of Figure 4 (Light Gray and Dark Gray). These distortions decrease when the distance between the two objects increases. This suggests that SDMs may not be the best LogOdds representation of discrete data for representing close objects. Of note, if linear operations are performed on the SDMs, the result will be a LogOdds but likely not an SDM, making its interpretation in terms of shape difficult. For Gaussian PAs on the other hand, the two distributions are not impacted by each other regardless of distance between the objects or the setting of the Gaussian filter. Thus, Gaussian PAs seem to better capture the shape of the two circles (see third row of Figure 4 (Light Gray and Dark Gray) as they were directly designed for the space of discrete probabilities.

The main advantage of the LogOdds framework is that if one wants to perform standard statistical analysis on shape, any form of PA can be used, as long as it is made up of valid discrete distributions.

3.2 Defining Rater-Specific LogOdds Maps

For many anatomical structures, manually tracing its boundary from an MR image is a serious challenge, as the MR signal does not provide enough contrast to clearly see the outline of the structure. This causes variations among expert segmentations of the same structure (see Figure 1). In this section, we show how PAs and thus LogOdds can be constructed to capture this variability.

Our representation makes use of the STAPLE algorithm presented in (Warfield et al., 2006). The method takes the set of binary maps traced by experts and turns them into SDMs. It then estimates a reference SDM via STAPLE based on the agreement between the experts' SDMs. We compute the performance parameters of each expert ε by calculating the mean, μ^(ε), and the variance, σ^(ε), of the voxel-wise difference between the SDM of the expert, 𝒟_(ε), and the reference SDM, 𝒟^(ℛ). The mean indicates an overall over- or under-estimation of the size of the structure by the expert, and the variance captures his ability to trace consistently. Given this model, the probability of a distance value ${\mathcal{D}}_x^{\left(\mathcal{E}\right)}$ at voxel x for expert ε is defined by the Gaussian distribution $\mathcal{N}({\mathcal{D}}_x^{\left(\mathcal{E}\right)}-{\mathcal{D}}_x^{\left(\mathcal{R}\right)}-{\mu }^{\left(\mathcal{E}\right)};{\sigma }^{\left(\mathcal{E}\right)})$.

Assuming that the voxels in the image are independently distributed, the probability of a voxel x being inside the object given 𝒟^(ε), μ^(ε) and σ^(ε) is then defined according to Bayes' rule as:

where $\Phi \left(y\right)\triangleq \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2}}\right)\right]$ is the Gaussian cumulative distribution and erf(·) the error function.

We can now interpret the map defined by the voxel entries $P({\mathcal{D}}_x^{\left(\mathcal{R}\right)}\ge 0\mid {\mathcal{D}}_x^{\left(\mathcal{E}\right)})$ as a rater-specific PA embedded in ${\mathbb{P}}_2^n$. A natural definition of the mapping function 𝒯 of a binary map B^(ε) drawn by an expert to LogOdds space would therefore be the LogOdds of the conditional probability:

Figure 5 shows the graph of the Gaussian cumulative functions computed from manual segmentations of the superior temporal gyrus by six experts (see Figure 1). The perfect segmentation would result in a step function (see light gray curve in Figure 5). All experts seem to perform similarly, except for expert D who is not only far from the truth but also unreliable, and expert E who shows great consistency and high accuracy. Similar observations can be made by looking at the corresponding LogOdds maps to the right of the graph, as the slope of the map indicates the overall accuracy of the expert. The map of expert E is close to a binary map indicating a high degree of agreement to the reference standard, whereas expert D shows a much softer LogOdds map.

3.3 Defining a time continuous atlas based on a finite number of samples

Neuroscientists often carry out longitudinal studies to better understand the aging process of a population. These studies are frequently defined by a set of subjects that have been scanned at different time points. We now explore the use of the LogOdds function for interpolating longitudinal data between time points.

This example is based on a longitudinal data set consisting of eight schizophrenic patients. Each patient has been scanned three times with an average separation of 14 months between the first and second scan, and an average 23 months separation between the second and third scan. For each time point, we generate a PA by first aligning the 24 MR images towards their central tendency using affine transformations computed from the population registration framework of (Zöllei et al., 2005). Afterward, we segment the gray matter from each aligned MR image using an atlas-based EM segmentation algorithm (Pohl et al., 2004). We then compute the PA of an anatomical structure at a given time point based on the overlap of the eight corresponding gray matter segmentations. The first row of Figure 7 shows a sample slice of the PA at three time points.

We can interpolate among atlases in either PA or LogOdds space. However, one should not perform linear operations directly on PAs, unless one restricts these operations to convex combinations (CC). As longitudinal data are generally composed of a few time points, only a very limited set of CCs can be applied to this type of data. Moreover, computing CC of PAs in the original space does not preserve the characteristic of certain Gaussian distributions over space as shown in Figure 6. In this example, the PA of a population is defined by a Gaussian distribution with mean A at time point 0 and B at time point 1. We therefore expect an interpolation between time point 0 and 1 to preserve the “hump” characterizing the distribution (see top graph, right column). This shape disappears at time point 0.5, when computing the CC of the two distributions within the PA space (see middle graph, right column). We can address this issue, however, by mapping the PAs into the LogOdds space and performing the CC there (see bottom graph, right column).

Most importantly, LogOdds give us the ability to use non-convex interpolation techniques. For longitudinal data, this provides us with a much richer class of interpolation functions than are available through CCs. An area of particular interest is the region around the thalamus (Figure 7 (a)) where an increase in volume can be observed over time in the PA. The second row of Figure 7 shows a magnified version of the PA of the thalamus region (a), and the corresponding linear (b) and quadratic spline interpolations (c). In the graphs, the z-axis and the intensity symbolize the probability, the x-axis is the time axis, and the y-axis represents the row of voxels highlighted by the black line in Figure 7 (a). Unlike (b), the quadratic interpolation is differentiable over time, which could enable us to extract additional parameters from the data such as the rate of change in the aging process of the population.

This completes the discussion of three examples that show the advantages of this representation over existing technologies. We first related SDMs to LogOdds and compared them to an alternative that was based on smoothing by spatial Gaussians. We then captured the uncertainty of boundaries by combining multiple segmentations of the same image to one LogOdds map. The last example described a framework for increasing the temporal resolution of PAs by interpolating the atlases within the LogOdds space.

4.1 Generating a deformable atlas via PCA on ${\mathbb{L}}_{M-1}^n$

We generate the deformable atlas by first turning a set of k manual segmentations {B⁽¹⁾,…,B^(k} into SDMs {ℱ^(l),…,ℱ^(k)}. Note that the segmentations have been aligned to each other through registering the corresponding images via (Zöllei et al., 2005). At this stage, our representation is very similar to that of (Tsai et al, 2003; Leventon, 2000; Yang et al., 2004; Golland et al., 2005). However, we interpret these SDMs as LogOdds maps and can thus embed the representation within a vector space. As such, we can perform PCA on the training set {ℱ⁽¹⁾,…,ℓ^(k} (Gentle, 1998). The deformable atlas is now defined by the eigenvector or modes of variation matrix U and $\stackrel{-}{\mathcal{F}}\triangleq {({\stackrel{-}{\mathcal{F}}}_1^T\cdots {\stackrel{-}{\mathcal{F}}}_{M-1}^T)}^T$, where ℱ̄_a is the mean vector of LogOdds of label a.

The deformable atlas (ℱ̄, U) encodes shapes $\stackrel{-}{\mathcal{F}}\left(\theta \right)\in {\mathbb{L}}_{M-1}^n$ within that atlas space by the expansion coefficient θ with ℱ(θ) = ℱ̄ + U · θ. We refer to the LogOdds maps of a label a defined by θ as

where matrix U_a is the entry in U corresponding to structure a. We define $U_{a,i}\in {\mathbb{L}}_1^n$ as the i^th eigenvector of structure a.

4.2 A Shape-Based Segmenter for MR Images

We now make use of the shape atlas by integrating it into a voxel-based classifier for the segmentation of anatomical structures with weakly visible boundaries in MR images. The method is an extension of a class of voxel-based classification methods (Wells et al., 1996; Van Leemput et al., 1999; Kapur, 1999; Marroquin et al., 2003; Pohl et al., 2004). These classifiers simultaneously estimate the image inhomogeneity β and determine the underlying label map 𝒯 based on the observed MR image ℐ. We extend this model by incorporating the expansion coefficients θ of the PCA model into the estimation process. We provide a full derivation of the new algorithm in Appendix A.

Briefly, the resulting EM-implementation is defined by two steps. The Expectation-Step (E-Step) computes the weights for each structure a and voxel x. The weights are defined by the product between the label map probability P(𝒯_x(a) = e_a|θ′) conditioned on the shape parameter θ′ and the intensity probabilities $P({\mathcal{I}}_x\mid {\mathcal{T}}_x\left(a\right)=e_a,{\beta }_x^{\prime })$ conditioned on voxel x being assigned to label a with image inhomogeneity ${\beta }_x^{\prime }$:

The Maximization Step (M-Step) estimates the inhomogeneities β′ and shape θ′ based on the weights 𝒲_x. The estimates for the shape parameters are the solution to the following MAP estimation problems:

The solution of Equation (5) depends on the definition of the conditional probability P(𝒯_x = e_a|θ) that captures the relationship between the shape parameters θ and the label map 𝒯. We now present three different interpretations of P(𝒯_x = e_a|θ) based on the shape atlas described in Section 4.1.

The first probabilistic model of P(𝒯_x = e_a|θ) is based on the assumption that the deformable shape atlas represents a family of level set functions whose zero level sets are the object's boundary. This representation reduces ℱ(θ)_a to a binary map via the Heaviside function:

A natural definition of the conditional probability is therefore

with normalization factor Z = 1 + Σ_{a′=1,…,M−1} ℋ(ℱ(θ)_a′,x) and M being the background. In general, Equation (6) defines a distribution over space with a very steep slope along the boundary of the structure. This causes the EM implementation to overemphasize the shape atlas compared to the image data, when computing the weights in the E-Step (see Equation (4)).

The second probabilistic model interprets the result of the PCA as members of the LogOdds space ${\mathbb{L}}_{M-1}^n$. According to Section 2.1, the inverse of the multinomial logit of the LogOdds entry ℱ(θ)_a,x defines the probability that voxel x is assigned to label a so that

We define the normalization Z ≜ 1 + Σ_{a′=1,…,M−1}^{eℱ(θ)_a′,x}. The conditional probability is now characterized with a more gradual slope along the boundary of the structures, allowing for more flexibility in determining the contour of the object.

The third probabilistic model is similar to the previous one, but the shape parameter θ = 0 is fixed. Thus, the PA for this model is defined as the multinomial logit of the mean LogOdds function ℱ̄

with Z ≜ 1 + Σ_{a′=1,…,M−1} ^{eℱ̄_a′,x}. This model ignores the modes of variations of the atlas and is thus quite similar to static PAs already proposed by (Van Leemput et al., 1999; Pohl et al., 2006).

We constructed three different EM implementations that only differ in the mapping of the shape atlas to PAs that is captured by the definition of P(𝒯_x = e_a|θ). The first and third mapping (Equation (6) + (8)) were influenced by representations of shape variations commonly used in the literature and the second mapping (Equation (7)) embodied the contribution of this paper. In the next section, we measure the robustness of the three implementations to provide the reader with a meaningful comparison of our LogOdds representation to well-established techniques. We note, however, that many other definitions of P(𝒯_x = e_a|θ) are possible.

4.3 The Accuracy of Bayesian Classifier

We now investigate the impact of the three implementations on the accuracy of the segmentation algorithm. We base the implementation EM-ℋ on the level set representation captured by Equation (6), the implementation EM-𝒫 on our new shape representation as defined by Equation (7), and the implementation EM-N (see also (Pohl et al., 2004)) on the static PA as described by Equation (8). All three implementations segment the caudate nucleus and the thalamus in 20 test cases. These two structures are of special interest for evaluation as they are characterized by very blurry boundaries in MR images. They are also characterized by different types of shapes, the thalamus being very round and the caudate more elongated.

We determine the quality of the automatic segmentations by comparing them to manual segmentations using a measure of overlap, the Dice coefficient (Dice, 1945). The score resembles the overlap between manual and automatic segmentation with a higher score given to those automatic segmentations that have greater overlap to the manual ones. The graph in Figure 8 shows the average Dice measures and standard error for the three implementations. If we interpret the manual segmentations as the gold standard then the average Dice score represents the accuracy of the implementation. The standard error is a measure of reliability with a small error indicating low fluctuation in performance of the approach. The discrepancy in performance is especially striking between EM-ℋ and EM-𝒫. For both structures, EM-𝒫 achieves a significantly higher average score (thalamus: 88.4±0.8%, caudate: 84.9 ±0.8%; mean DICE score ± standard error) than EM-ℋ (thalamus: 85.3 ±1.2%, caudate: 74.3 ±1.6%) as the range of scores defined by the means and standard errors of each structure do not overlap between the two implementations. EM-𝒫 also performs much better than EM-N(thalamus: 87.3 ±1.2%, caudate: 82.7 ±1.2%) with a significantly higher score for the caudate and better average score and lower standard error for the thalamus. We also note that EM-N achieves a much higher average score than EM-ℋ.

We think there are several reasons why EM-𝒫 is more accurate than the other algorithms, although their main drawback is that they do not properly capture variations within the population of shapes. EM-N only incorporates a static model of the shape variations, and is thus restricted to impose a shape constraint that only models the mean of the population. For the thalamus, this limitation only slightly reduces the accuracy of the approach (in comparison to EM-𝒫) because the shape of the structure hardly changes across the healthy population. The opposite, however, is true for the caudate, where the accuracy of the approach is much lower than reported for EM-𝒫. The caudate's elongated shape wraps around the ventricles, which can greatly differ in size.

The Heaviside function induces a shape model that is too strong compared to the image model, even when modes of variations are taken into account. The intensity information is especially important for segmenting the caudate as the elongated shape is more easily determined by the clearly visible boundary to the neighboring ventricles. Since the approach largely ignores the intensity information, it produces less reliable results than EM-𝒫. We believe EM-𝒫 to be a more flexible approach as it captures more variability and uncertainty, and thus only imposes the shape model when boundaries are weakly visible in the MR image.

In summary, we presented a statistical framework for the segmentation of anatomical structures in MR images. The framework is guided by the low-dimensional PCA shape model of Section 4.1 as the shape representation is turned into PA. We derived three different implementations that only differed within the probabilistic model. We ran each implementation on 20 test cases segmenting the thalamus and caudate. The segmentation algorithm based on our new representation EM-𝒫 performs much better than EM-ℋ that is based on a level set representation, or EM-N that uses a “conventional” PA such as described by (Van Leemput et al., 1999; Pohl et al., 2006).