Yunlong He
email: yunlong.he@insa-lyon.fr

Wanzhen Jiao

Yunfeng Shi

Jian Lian

Bojun Zhao
email: zhaobojun314@gmail.com

Wei Zou

Yuemin Zhu

Yuanjie Zheng
email: zhengyuanjie@gmail.com

Z Wei

Segmenting Diabetic Retinopathy Lesions in Multispectral Images Using Low-Dimensional Spatial-Spectral Matrix Representation

Keywords: Image segmentation, multispectral retinal images, spatial and spectral information, supervised feature learning

Multispectral Imaging (MSI) provides a sequence of en-face fundus spectral slices and allows for the examination of structures and signatures throughout the thickness of retina to characterize diabetic retinopathy (DR) lesions comprehensively. Manual interpretation of MSI images is commonly conducted by qualitatively analyzing both the spatial and spectral properties of multiple spectral slices. Meanwhile, there exist few computerbased algorithms that can effectively exploit the spatial and spectral information of MSI images for the diagnosis of DR. We propose a new approach that can quantify the spatialspectral features of MSI retinal images for automatic DR lesion segmentation. It combines a generalized low-rank approximation of matrices (GLRAM) with a supervised regularization term (SRT) to generate low-dimensional spatial-spectral representations using the feature vectors in all spectral slices. Experimental results showed that the proposed approach is very effective for the segmentation of DR lesions in MSI images, which suggests it as an interesting tool for assisting ophthalmologists in diagnosing, analyzing, and managing DR lesions in MSI.

I. INTRODUCTION

M Ultispectral imaging (MSI) refers to a non-invasive optical technique that facilitates the fine representation of physical objects by using a series of spectral bands, enabling observers to acquire information far beyond human vision. This technique has seen widespread use in research areas of remote sensing [START_REF] Davis | Remote sensing: the quantitative approach[END_REF], [START_REF] Laliberte | Multispectral remote sensing from unmanned aircraft: Image processing workflows and applications for rangeland environments[END_REF], dermatology [START_REF] Bouchard | Ultrafast multispectral optical imaging of cortical oxygenation, blood flow, and intracellular calcium dynamics[END_REF], histopathology [START_REF] Boucheron | Utility of multispectral imaging for nuclear classification of routine clinical histopathology imagery[END_REF], dentistry [START_REF] Chung | Multispectral nearir reflectance and transillumination imaging of teeth[END_REF]. It has recently received increased attention in retinal imaging, e.g., the Retinal Health Assessment (RHA, Annidis Health System Corp.) [START_REF] Li | In vivo study of retinal transmission function in different sections of the choroidal structure using multispectral imaging[END_REF], [START_REF] Li | Retinal and choroidal oxygen saturation of the optic nerve head in open-angle glaucoma subjects by multispectral imaging[END_REF]. Unlike conventional retinal imaging techniques which are limited to capturing retinal features visible to the human eye, MSI allows clinicians to visually examine both the spatial and spectral characteristics of structures throughout the thickness of retina. Specifically, MSI employs an extensive range of carefully selected, discrete monochromatic LED-sourced wavelengths to create a sequence of en-face fundus spectral slices, as shown in Fig. 1. Each slice reveals both the spatial distribution and lightabsorbing properties of structures from the internal limiting membrane through to the choroid [START_REF] Denniss | Relationships between visual field sensitivity and spectral absorption properties of the neuroretinal rim in glaucoma by multispectral imaging[END_REF]. With the ability to collect spatial and spectral information simultaneously, MSI provides a new manner to comprehensively characterize subtle, deep or overlapping lesions of retina, and thus can help to diagnose ocular diseases earlier than traditional funduscopy [START_REF] Khoobehi | Hyperspectral imaging for measurement of oxygen saturation in the optic nerve head[END_REF], [START_REF] Li | Snapshot hyperspectral retinal imaging using compact spectral resolving detector array[END_REF].

One of the most common ocular complications is Diabetic retinopathy (DR), which is caused by prolonged hyperglycemia damaging the retina and may result in severe vision loss and irreversible blindness [START_REF] Fong | Retinopathy in diabetes[END_REF], [START_REF] Yau | Global prevalence and major risk factors of diabetic retinopathy[END_REF]. Numerous epidemiological studies have shown that early diagnosis and treatment of DR can protect against diabetes-related visual impairment [START_REF] Klein | The wisconsin epidemiologic study of diabetic retinopathy: Ii. prevalence and risk of diabetic retinopathy when age at diagnosis is less than 30 years[END_REF]- [START_REF]Early photocoagulation for diabetic retinopathy: Etdrs report number 9[END_REF]. Retinal imaging provides the only avenue to directly visualize and quantify retinal structures for detecting early lesions indicative of DR, including microaneurysms, hard exudates, hemorrhages, etc. [START_REF] Abràmoff | Retinal imaging and image analysis[END_REF]. Assessment of the presence, size and location of these lesions assists ophthalmologists in determining disease severity in the early stages for supplying therapeutic interventions on time. Over the past decades, advances in retinal imaging have provided significant improvements in the clinical care and management of patients with DR [START_REF] Lombardo | Adaptive optics technology for high-resolution retinal imaging[END_REF]- [START_REF] Cole | Contemporary retinal imaging techniques in diabetic retinopathy: a review[END_REF]. As an emerging retinal imaging modality, MSI has recently been studied extensively for facilitating early and accurate diagnosis due to its potential diagnostic strengths [START_REF] Xu | A light-emitting diode (led)based multispectral imaging system in evaluating retinal vein occlusion[END_REF]- [START_REF] Bartczak | Applicability of led-based light sources for diabetic retinopathy detection in retinal imaging[END_REF].

Currently, visual estimation of both the spatial and spectral features of MSI images, which is based on the educated opinion of ophthalmologists, remains the reference standard for MSI-based pathology diagnosis. However, this process is tedious, subjective, and highly time-consuming because it requires clinical doctors to assess intensity variations across different slices and compare these variation properties between different locations, while unfortunately, a MSI sequence typically consists of more than ten spectral slices. As practicable solutions, computer-based algorithms for automatic assessment have attracted increasing attention thanks to their ability to provide objective and accurate assessments and reduce labor costs [START_REF] Styles | Quantitative analysis of multi-spectral fundus images[END_REF]- [START_REF] Ly | Multispectral pattern recognition reveals a diversity of clinical signs in intermediate age-related macular degeneration[END_REF]. An important task for computer-based diagnosis is the automatic lesion segmentation which can help the ophthalmologists to rapidly and reasonably detect DR lesions in MSI images, and thus to focus more time on evaluating disease aggressiveness and treatment response.

DR lesion segmentation aims to find the exact location and area of suspected lesions on retinal image. It is a highly challenging task in the area of retinal image analysis due to several obstacles, such as the high variability of lesion shapes and sizes, the presence of easily-confused components (e.g., optic disk, vessels, etc.), and unpredictable artifacts from nonuniform illumination and involuntary eye movements during image acquisition [START_REF] Nie | 3d deep learning for multi-modal imaging-guided survival time prediction of brain tumor patients[END_REF].

To address the difficulties associated with DR lesion segmentation, researchers have developed a variety of methods over the past few years [START_REF] Mookiah | Computer-aided diagnosis of diabetic retinopathy: A review[END_REF]. Typical segmentation methods include Fuzzy C-Means (FCM) clustering algorithms to divide image pixels into diverse clusters [START_REF] Osareh | Automated identification of diabetic retinal exudates in digital colour images[END_REF]- [START_REF]Automatic exudate detection for diabetic retinopathy screening[END_REF], region growing methods to form different image regions based on some uniformity criteria such as gray level and colour [START_REF] Sinthanayothin | Automated detection of diabetic retinopathy on digital fundus images[END_REF]- [START_REF] Fleming | Automated detection of exudates for diabetic retinopathy screening[END_REF], and mathematical morphology operations performed by analyzing geometrical structures of certain retina components [START_REF] Sopharak | Automatic detection of diabetic retinopathy exudates from non-dilated retinal images using mathematical morphology methods[END_REF]- [START_REF] Welfer | A coarse-to-fine strategy for automatically detecting exudates in color eye fundus images[END_REF]. Other methods focus on combining two or more abovementioned methods [START_REF] Streeter | Microaneurysm detection in colour fundus images[END_REF], [START_REF] Niemeijer | Automatic detection of red lesions in digital color fundus photographs[END_REF], or try to classify lesion and non-lesion pixels by using the well-known machine learning algorithms [START_REF] Gardner | Automatic detection of diabetic retinopathy using an artificial neural network: a screening tool[END_REF], [START_REF] Zhang | A svm approach for detection of hemorrhages in background diabetic retinopathy[END_REF]. In most of the methods, an essential factor in achieving a successful segmentation is the discriminative ability of features to identify the pathological clusters, regions, or pixels [START_REF] Faust | Algorithms for the automated detection of diabetic retinopathy using digital fundus images: a review[END_REF]. The features used in these methods involve handcrafted features predefined by human researchers (e.g., intensity difference and gradient [START_REF] Tariq | Automated detection and grading of diabetic maculopathy in digital retinal images[END_REF], [START_REF] Akram | Detection and classification of retinal lesions for grading of diabetic retinopathy[END_REF], local binary patterns (LBP) [START_REF] Mookiah | Evolutionary algorithm based classifier parameter tuning for automatic diabetic retinopathy grading: A hybrid feature extraction approach[END_REF], etc.) and learning-based features acquired by learning latent but discriminative representations from raw image data via machine learning [START_REF] Niemeijer | Automated detection and differentiation of drusen, exudates, and cotton-wool spots in digital color fundus photographs for diabetic retinopathy diagnosis[END_REF]- [START_REF] Sopharak | Machine learning approach to automatic exudate detection in retinal images from diabetic patients[END_REF]. Recently, convolutional neural networks (CNNs), a well-known deep learning architecture, have also shown impressive performance in extracting features that are robust to DR lesion variations [START_REF] Gulshan | Development and validation of a deep learning algorithm for detection of diabetic retinopathy in retinal fundus photographs[END_REF], [START_REF] Lam | Retinal lesion detection with deep learning using image patches[END_REF].

However, most of the existing techniques are restricted to operating on a single retinal image, making them difficult to simultaneously extract spatial and spectral features across multiple spectral slices. A possible solution to this problem is to use a CNN architecture in which three-dimensional (3-D) convolutional kernels could be used to cover the entire MSI stack, but existing 3-D CNN techniques are mainly designed for commonly used images (e.g. MRI brain images [START_REF] Nie | 3d deep learning for multi-modal imaging-guided survival time prediction of brain tumor patients[END_REF], video images [START_REF] Ji | 3d convolutional neural networks for human action recognition[END_REF], [START_REF] Wang | Non-local neural networks[END_REF], or hyperspectral images [START_REF] Chen | Deep feature extraction and classification of hyperspectral images based on convolutional neural networks[END_REF]). In addition, training effective CNN models typically requires a huge amount of well labelled training samples [START_REF] Shen | Deep learning in medical image analysis[END_REF], and this is especially evident when applying 3-D CNNs to MSI images since the complex relationships in the data must be learned [START_REF] Paoletti | A new deep convolutional neural network for fast hyperspectral image classification[END_REF]. Unfortunately, obtaining such a huge amount of labelled samples is not trivial in our case due to the complexity of manual interpretation of multiple retinal slices and the limited number of well-experienced experts in clinical practice.

In this paper, we propose a learning-based approach by exploiting both the spatial and spectral features of MSI retinal images. It first consists in representing the feature vectors of each pixel in all spectral MSI slices by a two-dimensional (2-D) spatial-spectral matrix. Then, we formulate the feature learning task in a framework of generalized low-rank approximations of matrices (GLRAM) [START_REF] Ye | Generalized low rank approximations of matrices[END_REF], which aims to compute low-dimensional and compact representations of a series of spatial-spectral matrices. To achieve more discriminative features, a supervised learning model designed for regression tasks in [START_REF] Zhen | Descriptor learning via supervised manifold regularization for multioutput regression[END_REF] is adjusted to incorporate the supervision of labels into the learning process. Finally, the 2-D spatial-spectral matrix on the LBP features is constructed [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] to well capture the spatial characteristics of each spectral image.

The rest of this paper is structured as follows. In Section II, we describe the proposed approach in detail. The experiments and results are reported in Section III, followed by the discussion in Section IV and the conclusion in Section V.

II. SPATIAL-SPECTRAL REPRESENTATION-BASED SEGMENTATION APPROACH

The idea of the proposed approach is to learn representative features by jointly using the spatial and spectral information of MSI images for pixel classification. For each sampled pixel, we extract its LBP feature vectors from all spectral slices and concatenate them into a single matrix. Then, the obtained matrices are fed into the feature learning algorithm that is formulated as GLRAM framework with a SRT to learn new low-dimensional representations. Finally, the learned representations are vectorised and used in a SVM classifier to segment healthy and lesion pixels across MSI images.

A. Problem statement

Given a sequence of retinal MSI spectral slices, each slice having m pixels. Supposing that there are n pixels (where 1 ≤ n ≤ m) randomly sampled from the lesion and background regions according to the manually delineated ground-truth, then we can extract n features {f 1 , f 2 , • • • , f n } as training samples, together with their corresponding labels

{l 1 , l 2 , • • • , l n }. For i = 1, 2, • • • , n, f i ∈ R d
denotes the ddimensional vector of features extracted from all spectral slices at pixel i, and the class label l i is usually a discrete variable (e.g., l i = 1 if pixel i is lesion and l i = -1 otherwise). Our goal is to generate compact and discriminative feature representation of each f i for pixel classification.

B. Matrix-based spatial-spectral representation

To combine the spectral and spatial information efficiently, we consider more natural but powerful matrix representations instead of using vectorized training samples. Assume that d = s × k, where s is the number of spectral bands and k is the number of features extracted from each band. Then for each pixel i, the features can be arranged to form a new 2-D matrix as follow:

M i = [f i,1 , f i,2 , • • • , f i,s ] ∈ R k×s , (1)
where M i is a spatial-spectral matrix representation instead of the training sample f i at pixel i. Fig. 2 illustrates the detailed implementation of M i adopted in our approach. As illustrated, features on each channel of MSI are first extracted by using the LBP operators [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. Then, LBP feature vectors are generated by vectorizing the local LBP image patches of size p (centred at pixel i) of LBP feature images from all channels. The LBP feature vectors are finally arranged into a single matrix to form M i . Note that for each pixel i, M i simultaneously contains its spatial characteristics in columns and its spectral properties in rows. The use of LBP is inspired by the previous finding that retinal images usually have varieties of self-similar patterns at different scales which can be considered texture signatures [START_REF] Mookiah | Evolutionary algorithm based classifier parameter tuning for automatic diabetic retinopathy grading: A hybrid feature extraction approach[END_REF].

C. Generalized low-rank approximation of matrices

In order to achieve compact and low-dimensional representations for these spatial-spectral matrices {M i } n i=1 , we formulate the feature learning task into the GLRAM [START_REF] Ye | Generalized low rank approximations of matrices[END_REF] framework, since the GLRAM can operate on a collection of matrix representations to reduce noise and compute accurate low-rank matrices, and moreover, it allows us to explore distinctive physical meanings residing in rows and columns of matrices to find optimal representations. Specifically, for each M i , the GLRAM aims to extract a low-rank approximation by solving the following minimization problem:

arg min U,V, {Ai} n i=1 1 n n i=1 M i -U A i V T 2 F , (2)
where U ∈ R k×r1 and V ∈ R s×r2 (where r 1 k and r 2 s) are two transform matrices to distinctively project a spatialspectral matrix M i into a new low-dimensional space, and A i denotes a coefficient matrix which is a low-dimensional representation of M i . Two constraints U T U = I r1 and V T V = I r2 are added with I r1 ∈ R r1×r1 and I r2 ∈ R r2×r2 to ensure that U and V have orthogonal columns. In our case, r 1 and r 2 are two pre-specified parameters which can be varied to explore the information respectively in rows and columns of M for determining the dimensions of A.

D. Incorporation of SRT into GLRAM

To achieve more discriminative representations of {A i } n i=1 , we include SRT into GLRAM to take advantage of prior knowledge from the label space. We first construct a symmetric similarity matrix W ∈ R n×n with the binary class labels [START_REF] He | Locality preserving projections[END_REF] with the parameter ϕ ∈ N + :

{l 1 , l 2 , • • • , l n }. The element W ij of W is computed using a heat kernel
W ij = e -l i -l j 2 ϕ , i, j = 1, 2, • • • , n. (3)
The diagonal elements of W are set to zeros, i.e., W ii = 0. In the present study, we set the parameter ϕ = 4 for class labels (either -1 or 1). With the weights W, we minimize i,j

A i -A j 2 F W ij . (4)
Note that for the binary class labels, the maximum weight W ij = 1 can be achieved if and only if l i = l j . This choice of W ij imposes a penalty if the class label l i is different from l j . In other words, minimizing Eq. ( 4) attempts to ensure that the learned low-rank approximations with the same class label (either lesion or healthy) tend to be spatially similar in terms of Frobenius norm. Thus, the discrimination of new low-rank representations can be improved. Combining Eq. ( 2) and Eq. ( 4), we obtain our objective function:

arg min U,V, {Ai} n i=1 1 n n i=1 M i -U A i V T 2 F + λ i,j A i -A j 2 F W ij , (5)
where the first term (GLRAM) is to seek a low-rank approximation A for each input matrix M , while the second term introduces the supervision of labels to improve the discriminative ability of A with λ ∈ (0, ∞) being a parameter for controlling the trade-off between approximation accuracy and discriminative ability. The goal is to find the optimal solutions of U , V and {A i } n i=1 and use the obtained U and V to project spatial-spectral matrices {M i } n i=1 onto lowdimensional space for DR lesion segmentation.

E. Alternate optimization of U and V

To jointly solve for U , V and {A i } n i=1 in Eq. ( 5), we employ an alternate optimization strategy [START_REF] Ye | Generalized low rank approximations of matrices[END_REF] instead of directly solving the objective function. Specifically, the first term in Eq. ( 5) can be rewritten as follows based on the property of trace of matrices:

1 n n n=1 M i -U A i V T 2 F = 1 n n i=1 tr((M i -U A i V T )(M i -U A i V T ) T ) = 1 n ( n i=1 tr(A i A T i ) + n i=1 tr(M i M T i ) -2 n i=1 tr(U A i V T M T i )), (6)
where the term n i=1 tr(A i A T i ) comes from the fact that both U and V are orthogonal, and for any two matrices, tr(AB) = tr(BA). Moreover, the term

n i=1 tr(M i M T i ) is a constant because {M i } n
i=1 is given. Hence, the minimization of Eq. ( 6) is equivalent to minimizing

1 n n i=1 tr(A i A T i ) -2 n i=1 tr(U A i V T M T i ) . ( 7
)
Setting the derivatives of Eq. ( 7) with respect to A i to 0, the minimum of Eq. ( 7) is achieved only if

A i = U T M i V , for i = 1, • • • , n.
This implies that for each i, A i is uniquely determined by U and V with 5) and changing the sign of the optimization problem, we get an equivalent maximization problem as follows:

A i = U T M i V . Inserting A i = U T M i V into Eq. (
arg max U,V 1 n n i=1 U T M i V 2 F -λ i,j U T (M i -M j )V 2 F W ij . (8)
The objective function of Eq. ( 8) has no closed-form solution.

To find the optimal solution of U and V, we rewrite Eq. ( 8) as:

arg max U,V 1 n tr n i=1 U T M i V V T M T i U -λ tr   i,j U T (M i -M j )V W ij V T (M i -M j )U   . (9)
For a given V , the objective function of Eq. ( 9) is equivalent to maximizing tr(U T Ψ u U ), where

Ψ u = 1 n n i=1 M i V V T M T i -λ i,j (M i -M j )V W ij V T (M i -M j ) T . (10)
Note that the maximum is achieved only if U ∈ R k×r1 consists of the r 1 eigenvectors of matrix Ψ u corresponding to the r 1 largest eigenvalues. Similarly, for a given U , Eq. ( 9) is equivalent to the maximization of tr(V T Ψ v V ), where

Ψ v = 1 n n i=1 M i U U T M T i -λ i,j (M i -M j )U W ij U T (M i -M j ) T . (11)
Algorithm 1 Algorithm for learning spatial-spectral features from MSI images Input: A collection of MSI spectral images {I τ } s τ =1 and their corresponding ground-truth image G containing the manually degraded area of DR lesions. 1: Sampling n pixels from the lesion and background regions randomly based on G; 2: for each: spectral image I τ do 3:

Construct the corresponding LBP feature image F τ by computing the LBP code for all the pixels in it. 4: end for 5: for each: sampled pixel i 6:

Construct the spatial-spectral matrix M i by using {F τ } s τ =1 and Eq. (1). 7: end for 8: Compute the similarity matrix W using the n corresponding labels {l i } n i=1 ; 9: Initialize V (0) = (I r2 , 0) T and set n ← 1 ; 10: repeat 11:

Calculate the matrix Ψ u by using Eq. ( 10); 12:

Compute the r 1 eigenvectors {α j } r1 j=1 of Ψ u corresponding to the r 1 largest eigenvalues; 13:

U (n) = [α 1 , α 2 , α 3 , • • • , α r1 ]; 14:
Calculate the matrix Ψ v by using Eq. ( 11); 15:

Compute the r 2 eigenvectors {β j } r2 j=1 of Ψ v corresponding to the r 2 largest eigenvalues; 16:

V (n) = [β 1 , β 2 , β 3 , • • • , β r2 ]; 17:
n ← n + 1; 18: until The objective functions Eq. ( 9) converges Output: The projection matrices U and V .

The solution of V ∈ R s×r2 consists of the r 2 eigenvectors of matrix Ψ v corresponding to the r 2 largest eigenvalues.

The above observations provide an iterative strategy for obtaining the optimal solutions of U and V : optimizing U by fixing V , and optimizing V by fixing U . Specifically, for a given V , we obtain U by computing the eigenvectors of the matrix Ψ u . With the obtained U , V can be updated by computing the eigenvectors of the matrix v . We repeated this procedure until convergence. In our method, we employ singular value decomposition (SVD) [START_REF] Yang | Two-dimensional pca: a new approach to appearance-based face representation and recognition[END_REF] to solve the the standard eigen decomposition because the truncated SVD can achieve the best approximation of given matrices with the Frobenius norm [START_REF] Ye | Generalized low rank approximations of matrices[END_REF], [START_REF] Zhang | Low-rank matrix approximation with manifold regularization[END_REF]. Pseudo-code for the details of the proposed learning algorithm is given in Algorithm 1.

F. Implementation of lesion segmentation

Based on the obtained U and V , we adopt a support vector machine (SVM) classifier [START_REF] Friedman | The elements of statistical learning[END_REF] to perform training and testing for lesion segmentation. In the training process, we use the computed U and V to project the matrix representations {M i } n i=1 , M i ∈ R k×s onto the low-rank approximations {A i } n i=1 , A i ∈ R r1×r2 . These {A i } n i=1 are then vectorized as a series of (r 1 × r 2 )-dimensional feature vectors to train a SVM classifier with a Gaussian kernel. In the testing process, all pixels from the evaluation dataset are first arranged into

III. EXPERIMENTS AND RESULTS

A. Databases

The database for evaluating the proposed approach is composed of 50 sequences of MSI images acquired using an Annidis RHA T M instrument (Annidis Health Systems Corp., Ottawa, Canada), including 40 unhealthy sequences and 10 healthy sequences. These images are of oculus dexter (OD) or oculus sinister (OS) from 20 patients diagnosed with DR and 5 healthy subjects. They are in the format of DICOM with a bit depth of 16, and each image is of size 2048 × 2048. Each sequence comprises 11 spectral slices captured using wavelengths of 550, 580, 590, 620, 660, 690, 740, 760, 780, 810, and 850 nm, respectively. In each sequence, all the slices were registered using a feature-point-matching based approach [START_REF] Zheng | Joint alignment of multispectral images via semidefinite programming[END_REF] to eliminate the spatial misalignments introduced during the imaging process. From the unhealthy images, several types of DR lesions such as microaneurysm, retinal hemorrhages, hard exudates, cotton wool spots, or macular edema are documented. Two ophthalmologists manually drew the lesion areas by carefully observing and comparing different spectral slices. They were invited to mark all lesion pixels across multiple bands with the Paint tool in Microsoft Windows. The final manually annotated results were served as ground-truths in the experiments. In particular, the annotations from one of the two ophthalmologists were only used to validate the segmentation stability of our approach.

As color fundus photographs are now routinely used for clinical diagnosis of DR, the well known MESSIDOR database [START_REF] Decencière | Feedback on a publicly distributed image database: the messidor database[END_REF] constructed for evaluating computer-assisted diagnoses of DR was also used to assess the proposed approach. It contains 1200 color fundus images which were captured using a TopCon TRC NW6 Non-Mydriatic fundus camera with a 45 degree field of view. These images are packaged in three sets of 400 images, whose sizes are respectively 1440 × 960, 2240 × 1488, and 2304 × 1536 with 8 bits per color plane. For each image, two diagnoses were provided by medical experts: DR grade and risk of macular edema. We used only the DR grade in the present study. Table I provides the description of the DR grades with the number of images for each grade. To evaluate our approach, 128 images were randomly selected from the database (28 healthy images from grade 0 and 100 unhealthy images from grades 1-3) and then were resized to 800 × 600. From the unhealthy images, hemorrhages and microaneurysms regions were manually marked by the two ophthalmologists. These marks are taken as ground-truths.

B. Evaluation Metrics

Four quantitative metrics were employed to evaluate the segmentation performance in terms of sensitivity (Se), specificity (Spe), accuracy (Acc), and the area under the curve (AUC) [START_REF] Odstrcilik | Retinal segmentation by improved matched filtering: evaluation on a new high-resolution fundus image database[END_REF]. They were computed by comparing each resulting image with the corresponding ground-truth. Se and Spe refer to measures of effectiveness in correctly identifying lesion pixels and non-lesion pixels, respectively. Acc and AUC reflect the ability to correctly identify the total number of pixels. To be intuitive, the receiver operating characteristic (ROC) curve was also used for performance evaluation. It is a 2-D plot that allows us to analyze the trade-off between Se and Spe. A ROC curve closer to the top left corner means a better performance.

C. Experimental Setting

In our experiments, there are several parameters to control: λ, n itr (number of iterations), L s and L r (sampling points and radius of the LBP operator), p (size of local LBP image patch, as mentioned in Section II-B), r 1 and r 2 (reduced dimensionality). The parameter λ was empirically set to 0.8 for a trade-off between reconstruction accuracy and discriminative ability. We repeated our algorithm for a number of times, each time performing 15 -20 iterations. We found that the iterative solution converges fast within few iterations, and the experiments also show that the accuracy tends to be stable L s = 8 with a smaller L r . We thus simply chose n itr = 8, (L s , L r ) = (8, 2) to obtain satisfactory results in all experiments. In addition, the performance of the proposed approach highly relies on the choice of reduced dimensionality r 1 , r 2 and the patch size p, which will be investigated in Section III-D1.

For all experiments, a leave-one-out-cross-validation was performed. For the MSI database, we did the training and testing experiments 50 times, with each time leaving out one of the 50 sequences for training, and using only the omitted one sequence for testing. The final result is the average from the 50 runs. For the MESSIDOR database, the 128 color fundus images were divided into 8 sets. Then leave-one-out tests were performed on the 8 sets of images: at each time, 7 sets of images were used to train a model and the remaining one set was used for testing, and this step was repeated until all the 8 sets of images had been used for testing. All the experiments were performed on a single Linux machine with 3.6GHz CPU and 48GB of system memory (RAM).

D. Results

1) Effect of patch size and reduced dimensionality:

To illustrate the effect of patch size p and reduced dimensionality r 1 , r 2 on MSI-based DR lesion segmentation, we give in Fig. 3 the segmentation accuracies of the proposed approach with patch size p = 9, 17, 25, and 33 under various reduced dimensionalities. The parameters r 1 and r 2 were tested by varying one while keeping the other fixed, e.g., varying 1 by keeping r 2 = 5. The results show that the reduced dimensionality and the patch size have clear influence on segmentation performance. In Fig. 3 we can observe that the accuracy of segmentation was low when p = 9 or p = 33, and among the best when p = 17. Moreover, Figs. 3(a) and(b) also indicate that the proposed approach achieves the best performance when r 1 = 7 and r 2 = 5, respectively. We will use this optimal parameter setting in all experiments.

2) Joint spatial-spectral representation: The performance of our joint spatial-spectral representation on segmentation was evaluated by varying the number of spectral slices within each MSI sequence. Fig. 4 shows the ROC curves and the AUC values of the proposed approach on the MSI database with four different numbers of spectral slices, including 1 slice at 550 nm, 4 slices ranging from 550 nm to 620 nm, 8 slices ranging from 550 nm to 760 nm, and 11 slices ranging from 550 to 850 nm, respectively. The parameter r 2 used in this experiment was specially set to 1 to keep a uniform dimensionality of the low-rank approximations for comparison. As we can see from Fig. 4, the proposed approach achieves an AUC value of 0.819 on one spectral slice (550 nm). By increasing the number of spectral slices, the segmentation results are improved gradually. The best ROC curve and AUC value are obtained when all Fig. 4. ROC curves of the proposed approach when varying the number of spectral slices in a MSI sequence (550 nm, from 550 to 620 nm, from 550 to 760 nm, from 550 to 850 nm, respectively). the 11 spectral slices within the MSI sequences are used. That demonstrates the gradually increasing discriminative ability obtained by adding the spatial and spectral information jointly for feature representations.

3) Dimensionality reduction: In this experiment, the performance of the proposed approach in terms of dimensionality reduction was evaluated. Three different techniques for low-dimensional representations were performed on the MSI database, including SVD, GLRAM and the proposed approach, with segmentation results given in Fig. 5. The GLRAM algorithm was carried out by running our algorithm with the parameter λ = 0 to only keep the GLRAM term. For suitable comparison, the SVD algorithm was with the dimension reduction parameter d = r 1 × r 2 to maintain a uniform dimension in the low-dimensional space. In addition, the untransformed features ( Fig. 5) were obtained by vectoring all the matrix representations {M i } n i=1 of pixels as (k × s)dimensional feature vectors directly to train the SVM classifier. As we can observe from Fig. 5, the low-dimensional representations learned by SVD, GLRAM, and our approach produce better results than the untransformed features without dimensionality reduction. Moreover, the AUC values of the GLRAM and our approach are respectively 0.086 and 0.131 higher than that of the SVD algorithm. This indicates the high performance of dimensionality reduction obtained by applying two-sided transformations on matrices.

4) Incorporation of SRT into GLRAM: The performance of incorporating SRT into GLRAM was assessed by comparing our algorithm with the GLRAM as well as other techniques. From Fig. 5, we observe that the ROC curve of the proposed approach is above that of the GLRAM, and the AUC value of the proposed approach is higher than the GLRAM by 0.045. Table II summarizes the quantitative results of the untransformed features, SVD, GLRAM and the proposed approach. The proposed approach achieves the highest values of Se, Spe and Acc among all the methods. Fig. 6 illustrates the qualitative segmentation results of the proposed approach, the GLRAM, and the manual annotations on a MSI sequence. Note that the results of the proposed approach (Fig. 6(d)) are closer to the manually delineated ground-truth compared to the results of GLRAM (Fig. 6(c)).

5) Segmentation stability:

The segmentation stability of the proposed approach was tested by running our algorithm with different manual annotations. Fig. 7 shows segmentation results obtained by taking the manual annotations of the two ophthalmologists as ground-truth. The ROC curves corresponding to the first and the second ophthalmologists are plotted in red and blue dashed, respectively. Leave-one-out experiments were performed on the MSI database. In Fig. 7 we observe that the ROC curve for the first ophthalmologist is slightly but not significantly above the curve for the second ophthalmologist. Furthermore, the proposed approach obtains an AUC of 0.965 for the first ophthalmologist, and an AUC of 0.929 for the second ophthalmologist, both of which indicate high segmentation accuracy of our approach.

6) Lesion segmentation on color fundus images: The robustness of the proposed approach was evaluated on the MES-SIDOR database. In this experiment, the dimension parameter r 2 of our algorithm was set to 1, and the green channel of each fundus image was used as input (only the green channel is routinely adopted in the research community [START_REF] Azzopardi | Trainable cosfire filters for vessel delineation with application to retinal images[END_REF]). The well-known dimensionality reduction technique, principal component analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF], was also implemented for comparison. Table III gives the quantitative results of PCA, GLRAM, and our approach on the MESSIDOR database when taking different manual annotations as ground-truth.

results for the first ophthalmologist's annotations show that the proposed approach achieves the best segmentation performance in terms of Se, Acc, and AUC among all the methods, and with Spe = 0.914, which is only 0.02 lower than the highest score (0.916) achieved by GLRAM. For the second ophthalmologist's annotations, the Se, Spe, Acc, and AUC of the proposed approach are among the highest.

DISCUSSION

The proposed approach possesses several advantages in handling the problem of MSI-based DR lesions segmentation. First, the proposed 2-D spatial-spectral matrix provides a flexible yet solution to make use of both the spatial spectral information of multiple MSI images. Second, the GLRAM employs two-sided transformations to transform a collection of 2-D matrices into low-dimensional space, which allows us to exploit distinctive physical quantities (e.g., spatial characteristics and spectral properties) residing in rows and columns of matrices to learn compact descriptors. Third, the SRT encodes the intrinsic local geometrical structure of label space, enabling the learned patterns closely related to their class labels for generating more discriminative representations. Moreover, the LBP-based feature extraction and the regularization strategy embed prior spatial and label information into the learning process, making the proposed approach capable of learning effective representations with limited number of training samples.

In the experiments, the optimal values of the parameters p, r 1 and r 2 were used. The choice of p depends on several factors, such as image resolution, magnification and the area covered by DR lesions. A small value of p often leads to limited spatial information whereas increasing the value of p too much causes inaccurate information to be extracted from patches. The parameters r 1 and r 2 control the number of eigenvalues of Ψ u and Ψ v to form the models U and V , and thus influence the dimension of low-rank approximations A . The optimal value of r 1 depends on the inherent spatial characteristics of DR lesions on each spectral slice, and that of r 2 on the inherent spectral properties of DR lesions across different slices. A small value of r 1 or r 2 may lack inherent spatial or spectral information of DR lesions. Besides, according to Eqs. [START_REF] Li | Snapshot hyperspectral retinal imaging using compact spectral resolving detector array[END_REF] and [START_REF] Fong | Retinopathy in diabetes[END_REF], Ψ u and Ψ v are calculated from the spatial-spectral matrix M i which inevitably contains the information of easily-confused components (e.g., optic disk, vessels, etc.) or noise. Therefore, a large value of r 1 or r 2 may introduce eigenvectors dominated by other components or noise. By changing the value of p, r 1 and r 2 , our approach is flexible to be applied to various images and applications.

Another interesting point would be to compare the proposed approach with CNNs-based approaches. The power of CNNs lies in their deep neural networks that are shown to be pertinent for representing training samples and are able to find complex feature patterns from data. Moreover, convolutional kernels with different patch sizes used in 3-D CNN architectures are especially suitable for MSI data. In this direction, appropriate regularization strategies may need to be developed to avoid overfitting problem caused by the high complexity of 3-D input data and the limited number of training samples available in practice.

V. CONCLUSION

Retinal MSI generates a sequence of retinal spectral slices that allow ophthalmologists to characterize DR lesions comprehensively. We have proposed a new feature learning approach for DR lesions segmentation in MSI images. It is formulated as a generalized low-rank approximation of matrices with a SRT to learn low-dimensional spatial-spectral representations from the feature vectors in all spectral slices. Experimental results showed that the proposed approach is very effective for the segmentation of DR lesions in MSI images, which suggests it as an interesting tool for assisting ophthalmologists in diagnosing, analyzing, and managing DR lesions in MSI.

Fig. 1 .

 1 Fig. 1. A sequence of MSI images captured by Annidis RHA from a patient aged 40 and diagnosed with diabetic retinopathy. Ordered from left to right and from top to bottom, the first 11 images are captured with the wavelengths of 550, 580, 590, 620, 660, 690, 740, 760, 780, 810, and 850 nm, respectively.

Fig. 2 .

 2 Fig. 2. Implementation of LBP-based spatial-spectral matrix representation.

  : (N M = 0) and (N H = 0) 546 1 (0 < N M ≤ 5) and (N H = 0) 153 2 (5 < N M < 15) and (0 < N H < 5) 247 and (N V = 0) 3 (N M ≥ 15) or (N H ≥ 5) or (N V > 0) 254 N M , N H and N V denote the number of microaneurysms, hemorrhages and neovessels, respectively. the 2-D spatial-spectral matrices M and then transformed into low-rank approximations A by using U and V . The obtained A are vectorized as (r 1 × r 2 )-dimensional features and fed into the trained SVM model for classification.

Fig. 3 .

 3 Fig. 3. Segmentation performance with different patch sizes and dimensions. p is the patch size, r 1 the reduced dimensionality by U and r 2 the reduced dimensionality by V .

Fig. 5 .

 5 Fig. 5. ROC curves of the untransformed feature strategy, SVD, GLRAM and proposed scheme (GLRAM with SRT) in DR lesion segmentation on MSI images.

Fig. 6 .

 6 Fig. 6. Qualitative results of the proposed segmentation scheme with and without SRT. The lesion regions are indicated by the orange color added in the MSI spectral slice. The bottom row shows the close-ups cropped from the top images. From left to right: (a) Original spectral slice of 580 nm in Fig. 1. (b) Manual annotations (Ground-truth) from ophthalmologist. (c) Results of segmentation without SRT (GLRAM). (d) Results of segmentation with SRT.

Fig. 7 .

 7 Fig. 7. ROC curves of the proposed approach when taking respectively annotations provided by ophthalmologist # 1 and ophthalmologist # 2 as the ground-standard.

TABLE I RETINOPATHY

 I GRADE IN MESSIDOR DATABASE.

TABLE II QUANTITATIVE

 II RESULTS OF LESION SEGMENTATION WITH DIFFERENT STRATEGIES.

		Untransformed features	SVD	GLRAM	Proposed
	Se	0.610	0.696	0.857	0.925
	Spe	0.859	0.885	0.967	0.983
	Acc	0.853	0.880	0.962	0.981

TABLE III PERFORMANCE

 III OF LESION SEGMENTATION ON THE MESSIDOR DATABASE.

	Ground truth	Methods	Se	Spe	Acc	AUC
		PCA	0.792	0.874	0.870	0.840
	Opt. #1	GLRAM	0.778	0.916	0.909	0.850
		Proposed	0.849	0.914	0.911	0.885
		PCA	0.702	0.882	0.876	0.798
	Opt. #2	GLRAM	0.768	0.880	0.875	0.832
		Proposed	0.812	0.907	0.903	0.868

Opt. #1 and Opt. #2 denote the manual annotations provided by the first and the second ophthalmologists, respectively.

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. XX, NO. XX, MONTH YEAR

IEEE OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. XX, NO. XX, MONTH YEAR

ACKNOWLEDGMENTS

This work was made possible through support from the National Natural Science Foundation of China via grant 61572300, 81871508 and 61773246, Taishan Scholar Program of Shandong Province of China via grant TSHW201502038, Major Program of Shandong Province Natural Science Foundation via grant ZR2018ZB0419. DECLARATIONS OF INTEREST: NONE.