Measuring Spectral Inconsistency of Multispectral Images for Detection and Segmentation of Retinal Degenerative Changes

Multispectral imaging (MSI) creates a series of en-face fundus spectral sections by leveraging an extensive range of discrete monochromatic light sources and allows for an examination of the retina’s early morphologic changes that are not generally visible with traditional fundus imaging modalities. An Ophthalmologist’s interpretation of MSI images is commonly conducted by qualitatively analyzing the spectral consistency between degenerated areas and normal ones, which characterizes the image variation across different spectra. Unfortunately, an ophthalmologist’s interpretation is practically difficult considering the fact that human perception is limited to the RGB color space, while an MSI sequence contains typically more than ten spectra. In this paper, we propose a method for measuring the spectral inconsistency of MSI images without supervision, which yields quantitative information indicating the pathological property of the tissue. Specifically, we define mathematically the spectral consistency as an existence of a pixel-specific latent feature vector and a spectrum-specific projection matrix, which can be used to reconstruct the representative features of pixels. The spectral inconsistency is then measured using the number of latent feature vectors required to reconstruct the representative features in practice. Experimental results from 54 MSI sequences show that our spectral inconsistency measurement is potentially invaluable for MSI-based ocular disease diagnosis.

wavelengths. MSI differs from other retinal imaging instruments in that it exploits an extensive range of discrete monochromatic LED-sourced wavelengths (ranging from 550 nm to 950 nm) to create a collection of noninvasive en-face spectral sections, which highlight the anatomic and metabolic signatures throughout the thickness of retina and choroid. The long-wavelength light (beyond 600 nm) reveals properties of melanin while a slightly longer wavelength reflects liposfuscin 5 . In contrast, wavelenghths ranging from 533 nm to 850 nm are selectively absorbed by hemoglobin, melanin and macular pigments. An assessment of the spectral inconsistency allows ophthalmologist's early detection of pathologies. For example, early changes of retinal pigment epithelium (RPE) of RPE disruption are observed as normal structures in short wavelengths and a slight structural variation in pigment pattern in the red (620 nm, 660 nm, 690 nm and 740 nm) and infrared (760 nm, 780 nm, 810 nm and 850 nm) spectral slices 4 .
Visual estimation of the MSI spectral inconsistency remains the reference standard for MSI-based diagnostic pathology, with which the ophthalmologist qualitatively assesses image variations across different spectra and compares these variation properties between different locations. However, automated and quantitative assessment of the spectral inconsistency is increasingly required along with elevated usage of MSI. Automated tools are potentially capable of providing more prognostic and predictive markers, which might enable ophthalmologists to assess the aggressiveness of a disease in its early stage or review its response to therapy. In addition, automated analysis is not aimed at replacing but rather assisting the ophthalmologists, to not only increase diagnostic precision and inter-observer reliability but also reduce labor costs.
In this paper, we propose a method for measuring MSI spectral inconsistency based on an outlier detection framework, which can be used to detect retinal degenerations and segment the corresponding deteriorated regions. Specifically, we measure spectral consistency by extracting the common spectral properties of normal tissues and specify degenerations as outliers that bear inconsistent spectral properties with normal tissues. Mathematically, we define spectral consistency as the fact that the representative features of any pixel in any spectrum can be reconstructed by projecting linearly a unique pixel-specific latent feature vector with a spectrum-specific projection matrix. In contrast, the reconstruction of a spectrally inconsistent pixel requires more than one latent feature vector. The proposed method is founded on a probabilistic Gaussian mixture model and designed to find a MAP (maximum a posteriori) estimate of the projection matrix and the assignment to the latent feature vector(s) via a stochastic expectation-maximization (SEM) algorithm. One unique property of this algorithm lies in the fact that the latent feature vectors do not need to be explicitly resolved, leading to a robust and fast estimation.

Results
Our validation database is comprised of 54 MSI image sequences acquired by using an Annidis RHATM instrument (Annidis Health Systems Corp., Ottawa, Canada). These images are of oculus dexter (OD) or oculus sinister The spectral inconsistency measurement approach we present in this paper offers a score value for each pixel to indicate the probability (a value in [01]) of this pixel being degenerated. We carry out 3 experiments in order to validate the usefulness of the output of the proposed approach. First, we compare both qualitatively and quantitatively between the algorithm-measured spectral inconsistency and the manually delineated degenerated regions in the MSI sequence. Second, we show that the inclusion of spectral inconsistency as an additional feature can help boost the performance of the classification-based retinal degeneration segmentation algorithm. Third, we test the speed of and agreement between ophthalmologists in manual segmentation of retinal degenerations and compare between a case of using the spectral inconsistency as a guide and another case without using it.
To eliminate the spatial misalignment between MSI spectral slices, we employ the semidefinite programming-based joint alignment algorithm 8 to register the sequential images in each MSI sequence. This algorithm solves a low-rank and semidefinite matrix that stores all pairwise-image feature mappings by minimizing the total amount of point-to-point matching cost via a convex optimization of a semidefinite programming formulation. It takes in a complete consideration of the information aggregated by all point-matching costs and enables the entire set of pairwise-image feature mappings to be solved simultaneously and near-optimally. This algorithm has been shown to be extremely effective in aligning sequential multimodel images.

Agreement Between Spectral Inconsistency and Degenerated Retinal Regions.
To assess the agreement between the spectral inconsistency and the degenerated retinal regions, we randomly chose 34 MSI sequences (3 from healthy subjects and 31 from patients), and in each sequence from patients, two ophthalmologists manually drew the degenerated areas by carefully observing and comparing different spectral slices. We then ran our algorithm to measure the spectral inconsistency value for each pixel of each sequence. We then obtained a set of degenerative regions by establishing a threshold for the spectral inconsistency in a way that pixels with a value above this threshold are treated as degenerated and as healthy otherwise. By varying the threshold value, we plotted the receiver operating characteristic (ROC) curve using the ROCKIT algorithm 9 and calculated the area under the curve (AUC) values, as shown by the results in Fig. 2. In addition to the AUC value, we also found that the best segmentation (as shown in Fig. 3) accuracy of our algorithm is 0.73 and 0.75 when treating the manual delineations of ophthalmologist #1 and #2 as the gold standard, respectively. From Fig. 3, we can see that microaneurysm, retinal hemorrhages and hard exudates can be observed from this subject, and our algorithm detected them successfully.

Classification-based Degeneration Recognition by Using Spectral Inconsistency.
We trained a SVM binary classifier 10 with a Gaussian kernel to distinguish normal pixels from the degenerated ones. The classifier treats all spectral values of each pixel as its features. A leave-one-out cross-validation strategy is exploited, which involves treating each of the 34 MSI sequences as the validation set and the remaining as the training set and repeating on all these ways to cut the original set into validation and training. We then consider the spectral inconsistency measured by our algorithm as an additional feature of the SVM classifier in order to validate its value in this per-pixel classification. The AUC values obtained without/with the spectral inconsistency are 0.73/0.77 and 0.74/0.76 when treating ophthalmologist #1 (as shown in Fig. 4) and ophthalmologist #2 as the gold-standard, respectively. In addition, we found that the best accuracies without/with using the spectral inconsistency are 0.79/0.84 and 0.79/0.83 when treating ophthalmologist #1 and ophthalmologist #2 as the gold-standard, respectively.   Manual Retinal Degeneration Segmentation Guided by Spectral Inconsistency. The spectral inconsistency measured by our algorithm is also invaluable due to its guidance on ophthalmologist's diagnosis. This guidance is actually founded on the fact that the spectral inconsistency indicates the probability of the pixel being degenerated. By observing the spectral inconsistency map, ophthalmologists may locate and recognize the details of the degenerated area more easily and quickly.
To validate this benefit from the spectral inconsistency, the two ophthalmologists manually outlined the degeneration by observing all the MSI spectral slices and the generated spectral inconsistency map for the remaining 20 MSI sequences. The average time (averaged over these two ophthalmologists) spent on each sequence is reduced from approximately 25 minutes on the 34 sequences (for which the spectral inconsistency was not involved) to approximately 14 minutes on the 20 sequences (for which the spectral inconsistency was used). In addition, we found that the Dice's coefficient (averaged over these two ophthalmologists), which measures the overlapping ratio between the two ophthalmologists' segmentation, is improved from 0.84 (generated without the spectral inconsistency) to 0.91 (generated with the spectral inconsistency).

Discussion
MSI employs multiple wavelengths of light sources to capture images of the retina and choroid and provides the ophthalmologists with the identification, interpretation, diagnosis and management of ocular pathology via spectral dissections. MSI represents a major advance in ocular diagnostics 4 because the carefully selected spectral bands used by MSI are targeted to the clinically relevant structures and metabolic characteristics, particularly the ocular chromophores melanin and hemoglobin, as well as the fluorophore lipofuscin. MSI allows for an isolation, identification and interpretation of the early and subtle ocular pathologies that are often difficult when using a traditional instrument.
However, visual inspection-based diagnosis with MSI is absolutely not trivial although it remains the reference standard. In the process of visual inspection of the MSI images, the ophthalmologists observe and compare the spectral characteristics across different spectra in order to measure qualitatively the spectral inconsistency. This is rather difficult considering the fact that the number of spectral slices in an MSI sequence may be more than 10 (as in the data set we are using in this paper).
In this paper, we propose a generative model for quantitatively measuring the MSI spectral inconsistency, which is validated to be able to reveal retinal degenerations. This model is built on our proposed mathematical definition of spectral consistency, which considers a pixel to be spectrally consistent when its representative features at any spectrum can be reconstructed by a linear projection of a single latent feature vector with a spectrum-specific projection matrix. With this definition, spectral inconsistency means a requirement of more than one latent feature vectors in the linear reconstruction. One key advantage of the proposed model is that the latent feature vectors (which play an import role in measuring the spectral inconsistency) are free from being inferred explicitly, which makes our algorithm robust and fast.
We have several findings from the experimental results described in the last section. First, segmentation of retinal degenerations by establishing a threshold for the spectral inconsistency measured using our algorithm complies well with manual segmentation of experienced ophthalmologists. This also means that the spectral inconsistency is potentially an invaluable indicator of retinal degenerations. Second, the spectral inconsistency measurement can help boost the performance of a classification algorithm designed for distinguishing the degenerated pixels from the normal ones and therefore can be treated as a new biomarker, which is potentially useful for early detection of pathologies. Third, the spectral inconsistency measurement can also offer a guidance for the ophthalmologists on recognizing pathological areas from MSI images, which reduces human effort and increases human speed in analyzing the sequential slices of MSI.
Disregarding the concrete generative model (to be detailed in next section) for measuring the spectral inconsistency from MSI, we can explain the mystery of why the proposed spectral inconsistency measurement can help to recognize the pathological areas from MSI based on an application of a similar mathematical model with multi-view anomaly detection 11 . Anomaly detection means identifying outlier observations, which do not conform to an expected pattern in a dataset. Unsupervised anomaly detection assumes that the majorities of the observations in the provided dataset are normal and looks for observations that seem to least fit to the remainder (i.e., outliers) of the dataset. In our spectral inconsistency measurement based on the proposed generative model, we assume that most pixels in MSI are normal, and then focus our efforts on searching for the potential outlier (i.e., degenerated) pixels whose spectral properties are different from the normal ones. The resulted spectral inconsistency is therefore associated with the probability of being pathological. This unsupervised identification of retinal pathologies is potentially invaluable in early diagnosis and easier differentiation of occult or overlapping pathologies.

Methods
Our method is based on a basic assumption that normal retinal tissues are consistent over certain spectral properties, whereas retinal degenerations are inconsistent. We propose a generative model approach for extracting the consistent spectral properties of normal tissues and measuring the spectral inconsistency of each retinal pixel, with which we can detect and segment the regions bearing degenerative changes by establishing a threshold for the spectral inconsistency. Specifically, we define spectral consistency as the fact that the representative features at any spectrum of a pixel can be reconstructed by a linear projection of a single latent feature vector with a spectrum-specific projection matrix. In contrast, the reconstruction of spectrally inconsistent pixels requires more than one latent feature vectors. As a key advantage of the proposed model, the latent feature vector are not needed to be explicitly inferred, which makes the algorithm robust and fast.
SCIENTIfIC RepoRtS | 7: 11288 | DOI:10.1038/s41598-017-11730-y Probabilistic Model for Measuring Spectral Inconsistency. We are given a sequence of s retinal MSI spectral slices, each with n pixels. We then have a data matrix d denotes a vector of features extracted from all spectral slices at pixel i. Supposing that there are d j features extracted from the jth spectral slice, we then have = ∑ = d d j s j 1 . We use x ij to represent the vector of features extracted at pixel i and from the jth spectral band. Our goal is then to specify degenerative pixels that bear inconsistent observation features across different spectral slices compared with normal pixels.
We define mathematically that a pixel i is spectrally consistent if its representative features x ij at any spectrum j can be reconstructed by a linear projection of a single latent feature vector  ∈ z i l with a spectrum-specific projection matrix  ∈ × P j d l j , i.e., = x Pz ij j i . We can see from this definition that z i is shared in the reconstruction of all spectra for a spectrally consistent pixel i. At the same time, we assume P j is shared by all pixels. In contrast, for a pixel i that is spectrally inconsistent, different spectra might be reconstructed with different latent feature vectors, which are denoted by , where  + represents the set of positive integers. In the probabilistic model introduced in this paper, the probability of the data x ij for the ith pixel and jth spectrum is calculated with the following Gaussian mixture model where N(x|μ, Σ) denotes a Gaussian distribution with mean μ and covariance matrix Σ, are the mixture weights, σ is a precision parameter and I means an identity matrix. We generate the mixture weights w m with a Dirichlet process, which is accomplished by a stick-breaking process 12 with a concentration parameter γ. At the same time, we assume that σ and the latent feature vector z im are drawn from a Gamma distribution Gamma (α, β) and a Gaussian distribution N(0, I/(σr)), respectively.
To measure the spectral inconsistency of retinal MSI based on the probabilistic model in Eq (1), we employ the SEM algorithm 11 . This algorithm is an iterative method to find a MAP estimate of the projection matrix {P j } and a latent vector assignment where s ij indexes the latent vector in Z i used for reconstructing x ij . One strength of this algorithm lies in its property of being free from the requirement of an estimation of the latent feature vectors Z i , mixture weights W i and precision parameter σ because they can be marginalized out analytically.

E-
Step. The E-step of the SEM algorithm samples a latent vector assignment value for any pixel i and spectrum j while fixing the projection matrix P with the following probability function

M-
Step. The M-step of the SEM algorithm estimates the projection matrix P while fixing the latent feature assignment S with the following estimation equation, which is obtained by maximizing the logarithm of the joint likelihood of X and S given P, α, β, r and γ (here we omit the deduction details for economy). The SEM algorithm iterates the E-step and the M-step, which samples s ij for each pixel i and spectrum j using Eq. (2) and estimates the projection matrix P j for each spectrum j.
Measuring Spectral Inconsistency. For each pixel i, we measure its spectral inconsistency ψ i by computing the average probability of the event M i > 1 over all iterations of the SEM algorithm with the following equation Specifying Representative Features and Setting Parameters. The representative features are extracted to describe the local image pattern around each pixel and specified as the image gradient orientation 13 rather than the image intensity. As shown in several works 14,15 , image gradient orientations can boost the performance of image classification, multi-output regression and subspace learning. Specifically, we compute the image gradients and the corresponding gradient orientation at each pixel of each spectral slice. We then treat the vector formed by concatenating all gradient orientations with a local window (in a size such as 7 × 7) as the representative feature vector of this pixel.
We terminated the algorithm after 400 iterations and empirically set the parameters γ = 1.2, α = 1 and β = 1. We specified the length l of the latent feature vector with a cross-validation process. Specifically, we compared the manually drawn degenerative regions of retina with the algorithm's results generated with different l values and then chose the smallest number of l with the best accuracy. Data availability statement. The datasets analysed during the current study are available from the corresponding author on reasonable request.