Toward a Severity Index for ROP: An Unsupervised Approach

A. Feature Extraction and Representation

We carry out principal spanning forest algorithm [14] to trace the vessels and construct the vasculature tree. Figure 1 displays tracing result on an example image from our dataset. Based on the observation that an image can contain both tortuous and straight vessels, an image is represented with two-component Gaussian Mixture Model (GMM) distribution of its image features [15]. This notion is demonstrated in Figure 1 where we compute tortuosity for each vessel segment and fit a two-component GMM to the extracted feature values. The figure shows the histogram of tortuosity features along with the vessels that are associated with each component. As can be seen tortuos and straight vessels are well distinguished with the probabilistic representation.

Based on our previous analysis [16] on the performance of various features with the aforementioned GMM representation, we extract acceleration, which is a point-based image feature that accounts tortuosity of a point on a vessel.

B. Distance between Images

Since each image is represented with the GMM distribution of its image features, we compute the distance between samples based on the distance between their distributions. Let pⁱ(f) denote the GMM distribution of the feature in image i with means ${\mu }_1^i$ and ${\mu }_2^i$, variances ${\sigma }_1^i$ and ${\sigma }_2^i$, component priors ${\pi }_1^i$ and ${\pi }_2^i$. Then, $p^i(f)={\pi }_1^i\mathcal{N}(f;{\mu }_1^i,{\sigma }_1^i)+{\pi }_2^i\mathcal{N}({\mu }_2^i,{\sigma }_2^i)$. The inner product between the distributions of image i and image j is < pⁱ, p^j >= ∫ pⁱ(f)p^j(f)df. After some derivations, this is simplified as

Using the inner product in equation (1) Euclidean distance satisfies D²(pⁱ, p^j) =< pⁱ, pⁱ > + < p^j, p^j > −2 < pⁱ, p^j >.

C. Dimensionality Reduction

We employ six different manifold learning techniques that are well-established and popular in the literature: principal component analysis (PCA), multidimensional scaling (MDS) [9], Isomap [10], local linear embedding (LLE) [11], Laplacian eigen-maps [12] and diffusion maps [13]. These methods uses a distance matrix between samples and aims to preserve either local or global distance relations in the reduced dimension. We plug in our GMM-based distance matrix, that measures pairwise distance as the distance between the GMMs. We obtain 1-dimensional coordinates as a result of these techniques, which can be regarded as a severity index. Moreover, these coordinates can help to order images with respect to the amount of tortuosity observed in the vessels.

D. Evaluation

In order to compare the results of different manifold learning algorithms, we compute mutual information (MI), a measure of dependency between two random variables. Assume X and Y are two random variables. The MI between them is I(X; Y) = H(X) − H(X|Y) where H(·) and H(·|·) denote entropy and conditional entropy respectively.

We use the Kernel Density Estimation plug-in MI estimator. Assume there are N samples, x_i, i = {1, 2, ···, N}. With a Gaussian kernel, the probability density function is $p(x)=N^{-1}{\sum }_{i=1}^N{G}_{{\sigma }^2}(x-x_i)$ where $G_{{\sigma }^2}(x_i)={(\sqrt{2\pi }{\sigma }_i)}^{-1/2}\mathit{exp}\{-(x_i^2)/(2{\sigma }_i^2)\}$ and ${\sigma }_i^2$ is the variable kernel covariance for the ith sample. The plug-in estimate of entropy is $H(X)=-N^{-1}{\sum }_{i=1}^{N}logp(x_i)$. The probability density given class label is p(x|c_k) = Σ_i∈Ck G_σ² (x − x_i) where C_k is the set of indices for the samples from class c_k. Then, the estimate of the conditional entropy can be written as $H(X\mid Y)=-{\sum }_{k=1}^{k}p(c_k)m_k^{-1}{\sum }_{i=1}^{m_k}logp(x_i\mid c_k)$ where m_k is the number of samples in class c_k. We compute kernel size by Silverman rule [17] where a fixed isotropic kernel is used for all samples.

To quantify the consistency between resulting coordinates and diagnosis, we compute the MI between resulting reduced coordinates of each algorithm and the class labels.

A. Dataset Preperation

We use a dataset of 104 wide-angle colored retinal images acquired during routine clinical care and independently graded by three experts as plus, pre-plus, or normal, and compared to the clinical diagnosis. Gold standard is obtained using methods previously described [18]. There are 19 plus, 31 pre-plus and 54 normal samples in the dataset. The tracings are obtained from manual segmentations in order to avoid the bias an automated segmentation algorithm might introduce.

B. Pairwise Distance Matrix

Figure 2 shows the color map visualization of the pairwise distance matrix. The first 19 rows display the plus images, while the following 31 rows are pre-plus and the rest are normal samples. As can be seen the block diagonals are better visible between normal and plus samples while pre-plus class seems ambiguous as some of its samples are closer to normal and the others resemble to plus samples. Hence, discrete classification on the dataset is a challenging task given the distance matrix and a continuous index might help to order samples according to their severity.

C. Performance Evaluation

We run all methods given the GMM-based distance matrix and obtain a 1-D output. However, Isomap, LLE and Laplacian eigenmaps require a local neighborhood parameter k, while diffusion map method necessitates input parameters σ and α associated to neighborhood and regularization respectively. For parameter selection, we use jackknife estimate. We vary values for these parameters and pick the best one for each method that maximizes jackknife estimate of MI between resulting coordinates and the class labels.

Table I reports the MI values computed for each method given the distance matrix and the selected parameter values. The coordinates Diffusion map generate has the highest consistency with expert’s decision. Figure 5 shows the resulting coordinates for the methods. For better visualization we also plot the rank orders such that the rank of each image is obtained by sorting the resulting 1-D coordinates. Class labels are indicated by different colors. As can be seen, for most of the methods pre-plus samples are mixed into the normal and plus classes, which is also consistent with the distance matrix visualization in Figure 2 where pre-plus samples have a diverse range of distance relation with other classes. Diffusion maps method provides a relatively better discrimination between classes. An interesting observation is that most of the methods tend to keep normal samples close to each other while plus and pre-plus samples are sparsely distributed. This might be due to higher intraclass connectivity of normal samples compared to the other classes. Moreover, this also shows that making a decision of pre-plus or plus is more challenging then determining a normal sample.

Figure 3 displays some sample images along with the ranking coordinates of diffusion maps algorithm. The numbers above each image show the number of expert decisions plus disease (+), pre-plus (±), and normal (−). We notice that from left to right, the amount of tortuous vessels increase. Furthermore, plus disease decisions among experts increase from left to right.

Figure 4 shows the severity ranking obtained with diffusion map versus that obtained with expert consensus (extracted from pairwise comparisons of images by five experts). While the unsupervised ranking exhibits some level of correlation with experts, we anticipate much better results when, in the future, we develop severity indices that consider labels in the form of diagnosis (class label) and relative severity (pairwise comparisons).