Combining Registration Errors and Supervoxel Classification for Unsupervised Brain Anomaly Detection

Abstract

Automatic detection of brain anomalies in MR images is challenging and complex due to intensity similarity between lesions and healthy tissues as well as the large variability in shape, size, and location among different anomalies. Even though discriminative models (supervised learning) are commonly used for this task, they require quite high-quality annotated training images, which are absent for most medical image analysis problems. Inspired by groupwise shape analysis, we adapt a recent fully unsupervised supervoxel-based approach (SAAD)—designed for abnormal asymmetry detection of the hemispheres—to detect brain anomalies from registration errors. Our method, called BADRESC, extracts supervoxels inside the right and left hemispheres, cerebellum, and brainstem, models registration errors for each supervoxel, and treats outliers as anomalies. Experimental results on MR-T1 brain images of stroke patients show that BADRESC outperforms a convolutional-autoencoder-based method and attains similar detection rates for hemispheric lesions in comparison to SAAD with substantially fewer false positives. It also presents promising detection scores for lesions in the cerebellum and brainstem.

The authors thank CNPq (303808/2018-7), FAPESP (2014/12236-1) for the financial support, and NVIDIA for supporting a graphics card.

A Appendix

Image Foresting Transform

The Image Foresting Transform (IFT) is a methodology for the design of image operators based on optimum connectivity [11]. For a given connectivity function and a graph derived from an image, the IFT algorithm minimizes (maximizes) a connectivity map to partition the graph into an optimum-path forest rooted at the minima (maxima) of the resulting connectivity map. The image operation resumes to a post-processing of the forest attributes, such as the root labels, optimum paths, and connectivity values. IFT has been successfully applied in different domains, such as image filtering [10], segmentation [6, 22, 36], superpixel segmentation [5, 24, 42], pattern classification [2, 29], and data clustering [27, 32]. This appendix presents preliminary concepts and introduces the IFT algorithm.

Preliminary Concepts

Image Graphs. A d-dimensional multi-band image is defined as the pair $\hat{I} = (D_I, \vec {I})$, where $D_{I} \subset Z^{d}$ is the image domain—i.e., a set of elements (pixels/voxels) in $Z^{d}$—and $\vec {I} : D_I \rightarrow \mathbb {R}^{c}$ is a mapping function that assigns a vector of c intensities $\vec {I}(p)$—one value for each band (channel) of the image—to each element $p \in D_I$. For example, for 2D RGB-color images: $d=2$, $c=3$; for 3D grayscale images (e.g., MR images): $d=3$, $c=1$. We represent a segmentation of $\hat{I}$ by a label image $\hat{L} = (D_I, L)$, wherein the function $L : D_I \rightarrow \{0, 1, \cdots , M\}$ maps every voxel of $\hat{I}$ to either the background (label 0) or one of the M objects of interest.

Most images, like the ones used in this paper, typically represent their intensity values by natural numbers instead of real numbers. More specifically, $\vec {I} : D_I \rightarrow [0,2^{b}-1]$, where b is the number of bits (pixel/voxel depth) used to encode an intensity value.

An image can be interpreted as a graph $G_I = (D_I, \mathcal{A})$, whose nodes are the voxels and the arcs are defined by an adjacency relation $\mathcal{A} \subset D_I \times D_I$, with $\mathcal{A}(p)$ being the adjacent set of a voxel p. A spherical adjacency relation of radius $\gamma \ge 1$ is given by

$$\begin{aligned} \mathcal{A}_{\gamma }: \{(p, q) \in D_I \times D_I, \left\| q - p \right\| \le \gamma \}. \end{aligned}$$

(3)

The image operators considered in this paper use two types of adjacency relations: $\mathcal{A}_1$ (6-neighborhood) and $\mathcal{A}_{\sqrt{3}}$ (26-neighborhood), as illustrated in Fig. 12.

Fig. 12.

Examples of adjacency relation for a given voxel p (red). (Color figure online)

Paths. For a given image graph $G_I = (D_I, \mathcal{A})$, a path $\pi _q$ with terminus q is a sequence of distinct nodes $\langle p_1, p_2, \cdots p_k \rangle$ with $\langle p_i, p_{i+1} \rangle \in \mathcal{A}$, $1 \le i \le k - 1$, and $p_k = q$. The path $\pi _q = \langle q \rangle$ is called trivial path. The concatenation of a path $\pi _p$ and an arc $\langle p, q \rangle$ is denoted by $\pi _p \cdot \langle p, q \rangle$.

Connectivity Function. A connectivity function (path-cost function) assigns a value $f(\pi _q)$ to any path $\pi _q$ in the image graph $G_I = (D_I, \mathcal{A})$. A path $\pi ^{*}_q$ ending at q is optimum if $f(\pi ^{*}_q) \le f(\tau _q)$ for every other path $\tau _q$. In other words, a path ending at q is optimum if no other path ending at q has lower cost.

Connectivity functions may be defined in different ways. In some cases, they do not guarantee the optimum cost mapping conditions [9], but, in turn, can produce effective object delineation [26]. A common example of connectivity function is $f_{max}$, defined by

$$\begin{aligned} \begin{aligned} f_{max}(\langle q \rangle )&= \left\{ \begin{array}{ll} 0 &{} \text {if }q\in \mathcal {S}, \\ +\infty &{} \text {otherwise.} \end{array} \right. \\ f_{max}(\pi _p\cdot \langle p,q \rangle )&= \max \{f_{max}(\pi _p), w(p, q)\}, \end{aligned} \end{aligned}$$

(4)

where w(p, q) is the arc weight of $\langle p, q \rangle$, usually estimated from $\hat{I}$, and $\mathcal {S}$ is the labeled seed set.

Fig. 13.

Multi-object image segmentation by IFT. (a) Axial slice of a brain image with seeds $\mathcal {S}_0$ for the background (orange), $\mathcal {S}_1$ for the right ventricle (red), and $\mathcal {S}_2$ for the left ventricle (green). (b) Gradient image for (a) that defines the arc weights for seed competition. Arcs have high weights on object boundaries. (c) Resulting segmentation mask for the given seeds and arc weights. Red and green voxels represent object voxels, whereas the remaining ones are background. (Color figure online)

1.1 The General IFT Algorithm

For multi-object image segmentation, IFT requires a labeled seed set $\mathcal {S} = \mathcal {S}_0 \cup \mathcal {S}_1 \cup \cdots \mathcal {S}_M$ with seeds for object i in each set $\mathcal {S}_i$ and background seeds in $\mathcal {S}_0$, as illustrated in Fig. 13. The algorithm then promotes an optimum seed competition so that each seed in $\mathcal {S}$ conquers its most closely connected voxels in the image domain. This competition considers a connectivity function f applied to any path $\pi _q$.

Defining $\varPi _q$ as the set of all possible paths with terminus q in the image graph, the IFT algorithm minimizes a path cost map

$$\begin{aligned} C(q) = \underset{\forall \pi _q \in \varPi _q}{min}\{f(\pi _q)\}, \end{aligned}$$

(5)

by partitioning the graph into an optimum-path forest P rooted at $\mathcal {S}$. That is, the algorithm assigns to q the path $\pi ^{*}_q$ of minimum cost, such that each object i is defined by the union between the seed voxels of $\mathcal {S}_i$ and the voxels of $D_I$ that are rooted in $\mathcal {S}_i$, i.e., conquered by such object seeds.

Algorithm 2 presents the general IFT approach. Lines 1–7 initialize maps, and insert seeds into the priority queue Q. The state map U indicates by $U(q) = White$ that the voxel q was never visited (never inserted into Q), by $U(q) = Gray$ that q has been visited and is still in Q, and by $U(q) = Black$ that q has been processed (removed from Q).

The main loop (Lines 8–20) performs the propagation process. First, it removes the voxel p that has minimum path cost in Q (Line 9). Ties are broken in Q using the first-in-first-out (FIFO) policy. The loop in Lines 11–20 then evaluates if a path with terminus p extended to its adjacent q is cheaper than the current path with terminus q and cost C(q) (Line 13). If that is the case, p is assigned as the predecessor of q, and the root of p is assigned to the root of q (Line 14), whereas the path cost and the label of q are updated (Line 15). If q is in Q, its position is updated; otherwise, q is inserted into Q. The algorithm returns the optimum-path forest (predecessor map), root map, path-cost map, and the label map (object delineation mask).