Brain-inspired algorithms for retinal image analysis

Retinal image analysis is a challenging problem due to the precise quantification required and the huge numbers of images produced in screening programs. This paper describes a series of innovative brain-inspired algorithms for automated retinal image analysis, recently developed for the RetinaCheck project, a large-scale screening program for diabetic retinopathy and other retinal diseases in Northeast China. The paper discusses the theory of orientation scores, inspired by cortical multi-orientation pinwheel structures, and presents applications for automated quality assessment, optic nerve head detection, crossing-preserving enhancement and segmentation of retinal vasculature, arterio-venous ratio, fractal dimension, and vessel tortuosity and bifurcations. Many of these algorithms outperform state-of-the-art techniques. The methods are currently validated in collaborating hospitals, with a rich accompanying base of metadata, to phenotype and validate the quantitative algorithms for optimal classification power.


Introduction
Diabetes is reaching epidemic proportions worldwide, especially in Asia due to fast lifestyle changes and genetic factors.Today 11.6 % of the Chinese population is estimated to have diabetes-2.Diabetic retinopathy is also the main cause of newly formed blindness in the working population, leading to high societal costs.Early detection is the key to prevention and successful treatment of these forms of blindness.However, many cases still go unnoticed and are not treated in time, especially in rural areas.
To set up a screening program for the detection of early signs of diabetic retinopathy (DR), glaucoma and age-related macula degeneration (AMD), a Sino-Dutch consortium was formed and the project RetinaCheck was defined with the following four phases: (1) development of innovative algorithms for automated and quantitative detection of relevant bio-markers, (2) set up a significant validation study correlating the imaging data with relevant clinical metadata, and (3) roll out a screening infrastructure in the province of Liaoning, Northeast China, and (4) make a sustainable commercial infrastructure.The partners include the departments of Biomedical Engineering at Eindhoven University of Technology (TU/e, the Netherlands) and Northeastern University (NEU, Shenyang, China), which develop the computer-aided diagnosis (CAD) software, the clinical partners He University Eye Hospital (HUEH, Shenyang, China) and China Medical University (Shengjing Hospital, Shenyang, China) (Fig. 1), and the fundus camera manufacturer i-Optics (the Hague, the Netherlands).
This paper reports on phase 1 only, focusing in detail on algorithm design.The very important validation phase is ongoing, and will be reported in a forthcoming paper.
Due to its low cost, wide availability, high resolution and ease of use, optical imaging of the retinal fundus is often the method of choice for screening applications.Today non-mydriatic cameras (no pupil dilation) give good quality and high-resolution images.Quantitative measurement and analysis of geometric attributes of the retinal vasculature is Extensive clinical metadata are co-recorded for validation highly informative for diagnosis and monitoring of diabetes at early stages [41].The many extracted imaging biomarkers in this project will also be exploited for the early detection and monitoring of other diseases, like glaucoma, age-related macula degeneration (AMD), and cardiovascular and neurodegenerative diseases.
Several other early DR signs can be measured, such as nerve damage in the cornea with confocal laser microscopy, or changes in retina neural tissue layer thickness with optical coherence tomography (OCT), but these methods are more costly and more labor intensive, especially given the projected huge-scale screening and the much lower availability of OCT in China.
This paper gives an overview of the current brain-inspired applications in the RetinaCheck project.The paper mainly focuses on vessel analysis, while a forthcoming paper will focus on the automated detection of micro-aneurysms and background diagnostic features, as drusen and exudates.
The paper is organized as follows: first a short introduction is given to the physiological evidence of multi-scale and multi-orientation processing in the visual system in Sect.2, and a theoretical model for brain-inspired multi-orientation computing in Sect.3. Then a series of recently developed computer-aided diagnosis (CAD) algorithms is discussed in Sect.4, mostly based on the mechanism of invertible multiorientation scores.The paper ends with a description of the current validation study in Sect. 5 and a conclusion in Sect.6.

Brain-inspired computer vision
Modern brain imaging techniques have revolutionized brain research: optical imaging methods such as voltage sensitive dye imaging, calcium intrinsic imaging and optogenetics, as well as structural and functional MRI techniques have revealed intricate details of brain function and structure, especially of the visual system, the best studied brain area today.Many models of the functional mechanisms in early vision have been proposed.This paper focuses on the geometric approach, starting from the multi-scale and multi-orientation structure found in the early stages of vision.

Multi-scale analysis
Multi-scale analysis is now well established [15,44,58].The center-surround receptive fields in the retina are the first multi-scale sampling step.Our model for V1's 'simple cells' (so called by Hubel and Wiesel) is that these have the function of multi-scale, regularized spatial differential operators, i.e., Gaussian derivatives, possibly up to fourth order [44,[58][59][60].The simultaneous sampling of the multi-scale differential structure leads to 'deep structure' analysis [43], with applications as edge focusing, hierarchical top points and SIFT and SURF keypoints.Nonlinear adaptive multi-scale analysis is developed in the rich field of geometry-driven diffusion [32,62], which introduced nonlinear PDEs and energy minimizing variational methods into this field of evolutionary geometric computing.

Multi-orientation analysis
Hubel and Wiesel [37] discovered that the receptive fields in cat's striate cortex have a strong orientation-selective property (see Fig. 2).A so-called cortical hyper-column, with the characteristic pinwheel structure of equi-orientation lines radiating from a central singularity, can be interpreted as a visual pixel computer, neatly decomposed into a complete set of orientations.
Synaptic physiological studies of horizontal pathways in cat striate cortex show that neurons with aligned receptive field sites excite each other on the cortex between different pinwheels, in an elongated area and over long distances, creating long-range contextual connections [2].Therefore, the visual system not only constructs the aforementioned score of local orientations, but also accounts for long-range context and alignment by excitation and inhibition a priori, which can be modeled by so-called association fields [15,20,27], and left-invariant PDEs and ODEs for contour enhancement and contour completion directly on the score [23,61].
This multi-orientation framework [14,15,21,31] also allows us to generically deal with crossings, as we will show with the application to vessel tracking [6,8,9,52] and segmentation [1,33,65].Moreover, due to the neat organization of image data on the Lie-group SE (2), we are able to design effective detection algorithms [6,7], geometric feature analysis techniques such as bifurcation detection/analysis [55] and Fig. 2 Left voltage sensitive dye measurement of the visual cortex of the tree shrew, exhibiting the columnar organization with the typical multi-orientation pinwheel structure.From [40].Right Orientation score model with orientation on the vertical axis, and cortical columns with horizontal contextual (association field) connections.From [66] local vessel curvature analysis [11], and enhancement methods [31,66].

Theory
Motivated by this organization, so-called orientation scores are constructed by lifting all elongated structures (in 2D images) along an extra orientation dimension, see Fig. 2. Similar to the perceptual organization of orientation in the visual cortex, a 2D orientation score is an object that maps each 2D position and orientation (x, y, θ) to a complex scalar.So the original 2D image domain can be extended to the score domain.A great advantage is that it can deal with multiple orientations per position, and the extra dimension enables new techniques for, e.g., contextual geometric reasoning and crossing-preserving enhancement [31].

The Euclidean motion group SE(2)
The domain R 2 S 1 (2 spatial dimensions, and one orientation dimension) of an orientation score can be identified with the Euclidean motion group SE(2), equipped with the group product for all g, g ∈ SE(2), and with R θ = cos θ − sin θ sin θ cos θ a counterclockwise rotation over angle θ .

Invertible orientation scores on SE(2)
An orientation score U f : SE(2) → C is obtained by convolving an input image f with a specially designed, anisotropic wavelet ψ: where R θ denotes a 2D counterclockwise rotation matrix and ψ θ (x) = ψ(R −1 θ x).Exact image reconstruction is achieved by where F[•] represents the unitary Fourier transform on L 2 (R 2 ) and M ψ is given by 2π 0 |F [ψ θ ] | 2 dθ .Well-posedness of the reconstruction is controlled by M [10,24].M : R 2 → R + is calculated by The function M provides a stability measure of the inverse transformation.Theoretically, reconstruction is well-posed, as long as 0 < δ < M < M <, where δ is arbitrarily small, since then the condition number of W is bounded by Mδ 1 (for a detailed discussion see [22]).
Wavelets that meet this criterion are the so-called cake wavelets described by [6,21,31].They uniformly cover the Fourier domain up to a radius of about the Nyquist frequency ρ n (see Fig. 3).They have the advantage over other oriented wavelets (such as Gabor wavelets) that they allow for a stable inverse transformation W * ψ from the orientation score back to the image.As such, cake wavelets ensure that no data-evidence is lost during the forward and backward transformation.The procedure for creating cake wavelets is illustrated in Fig. 3.In detail, the proposed wavelet takes the form with G s an isotropic Gaussian window in the spatial domain and (ρ, ϕ) T denote polar coordinates in the Fourier domain, i.e., ω = (ρ cos ϕ, ρ sin ϕ) T .The Fourier wavelet ψ(ω) = A(ϕ)B(ρ) has an angular component A(ϕ) resembling a 'piece of cake' (Fig. 3), and a Gaussian weighted radial component B(ρ): where B k (x) denotes the kth order B-spline, N θ is the number of samples in the orientation direction and s θ = 2π/N θ is the angular step size.The function B(ρ) is a Gaussian multiplied with the Taylor series of its inverse up to order N to enforce faster decay.The parameter t is given by From [6] with the inflection point ρ that determines the bending point of B(ρ).
As depicted in Fig. 3, the symmetric real part of the kernel picks up lines, whereas the asymmetric imaginary part responds to edges.In the context of vessel filtering, the real part is of primary interest, whereas vessel tracking in the orientation score additionally makes use of the imaginary part [6].Elongated structures involved in a crossing can now be disentangled for crossing-preserving analysis.

Left-invariant Gaussian derivatives in SE(2)
Theoretically, because of the curved geometry of orientation scores, it is wrong to take the derivatives in orientation scores using the {∂ x , ∂ y , ∂ θ } derivative frame (we use shorthand notation are used in SE (2).The ∂ ξ and ∂ η are the spatial derivatives tangent and orthogonal to the orientation θ .It is important to mention that not all the left-invariant derivatives commute, e.g., [22,31].Mathematically this is due to the fact that rotations and translations do not commute.We have the following commutators in the Lie-algebra of SE(2): Suitable combinations of derivatives have been widely used to pick up differential invariant features like edges, ridges, corners and so on [59].However, obtaining derivatives directly is an ill-posed problem.Therefore, we regularize the orientation scores via convolutions with Gaussian kernels the spatial Gaussian distribution G σ s : R 2 → R + must be isotropic to preserve commutator relations of the SE(2) group for scales σ s > 0, i.e., to preserve left-invariance.Detailed numerical approaches for linear left-invariant diffusions on SE(2) have been developed in [66].

Computer-aided diagnosis algorithms
The pipeline of functions in the developed CAD system is the following: Each of the steps will be discussed below in the respectively numbered sections.

Automatic quality assessment
Sometimes retinal images have a low quality, e.g., due to cataract or other pathology.To exclude non-diagnostic images automatically in the high-volume screening process, a method for image quality verification is developed [26], based on [50].The supervised method is based on the assumption that sufficient image structure according to a pre-defined distribution must be available.Geometric differential invariants (expressed in 2D gauge coordinates [58]) up to second order and at 4 different scales are used to develop response vectors.Different combinations of features are trained by different classifiers on 100 normal and 100 lowquality images.The ground truth for normal or low quality images was specified by two expert ophthalmologists.Combining the image structure clusters with RGB color histogram features, the Random Forest classifier proved to be the best classifier, with a performance of 0.984 area under the curve (AUC) of the receiver operator characteristic (ROC), with 0.91 accuracy rate.

Masking and normalization
All pixels outside the camera field-of-view are masked, according to the camera type.Retinal images often suffer from non-uniform illumination and varying contrast, which may affect the later detection process.We exploit the luminosity and contrast normalization procedures proposed by Foracchia et al. [28].

Optic nerve head detection
The optic nerve head (ONH) or optical disk is a key landmark in retinal images, and many methods for detection have been proposed [4,18,45,47,54,64].For an extensive and recent overview see [51].Correct segmentation of the ONH and its rim is an important biomarker for glaucoma identification [38]), and for establishing a metric for regions-of-interest on the retina [36].
Conventional fundus images show the ONH as a bright disk-like feature.However, scanning laser ophthalmoscopy (SLO) cameras generally show dark regions with considerably less contrast.In the presence of large pathologies, classical approaches typically show decreased performance, or fail.Better performance is obtained by including contextual information, e.g., by incorporating the characteristic pattern of large blood vessel arches in the upper and lower retina emerging from the ONH [45,46,64], but this comes at higher computational costs.
Exploiting the extra dimensions of orientation scores, we have developed a template matching application via crosscorrelation including local orientations [7,10] (see also [16]).We can now match patterns of orientation distributions, rather than pixel intensities.
The templates are constructed from (a) a disk filter, (b) a model of the vessels radiating from the ONH, and (c) through minimization of a suitable energy functional.In Ref. [7] we found that, in contrast to model-based templates [10] (e.g., disk shapes or vessel indicator functions), the best performing templates are trained by the minimization of an energy functional, consisting of a term where the desired inner product responses are trained, and a regularization term that ensures stability and smoothness of the template.For 2D spatial image templates t, we minimize where fi is one of P (P = 100) normalized positive patches with ONH (with y i = 1), or negative patches without ONH (with y i = 0) with typically λ = 10 −1.5 .The regularization term enforces smoothness by punishing the squared gradient magnitude ∇t 2 , and prevents over-fitting.Examples of both positive and negative training patches are given in Fig. 4.
For the orientation score template T , a similar functional is minimized: with typically λ = 10 and D θθ = 10 −3.5 and with the left- . The templates are represented in a B-spline basis, allowing for efficient and accurate optimization of the energy functionals.The B-spline coefficients are solved by a conjugate gradient approach.
The orientation score SE(2) template method outperforms state-of-the-art methods on publicly available benchmark databases, as it correctly identifies the ONH in 99.7 % of 1737 images of the well-known MESSIDOR, STARE (with a wide variety of pathological images) and DRIVE databases.For more details, see [7,10].

Vessel enhancement filter
Vessels can be enhanced by oriented filters [5,29], e.g., constructed from second-order Gaussian derivative operators at appropriate scales and orientations.To solve segmentation problems at complex structures like crossings, we have designed new filters in the orientation score matching the vessel profile as second-order left-invariant Gaussian differential operators perpendicular to the corresponding orientation [65], as described in Sect.3.
With Gaussian filters, the maximum response occurs at σ = r/ √ 2, where r represents the radius of the vessel caliber [44].Typically vessel calibers of the DRIVE and STARE databases range from 2 to 14 pixels, so we sample the spatial scales σ s as S = {0.7,1.0, 1.5, 2.0, 2.5, 3.5, 4.5} and angular scale σ o = π/5.We exploit N o = 36 orientations between 0 and π .Figure 5 shows our segmentation results on the DRIVE and STARE databases.From Fig. 6 we can see that the proposed orientation score based multi-scale filters show much better structure preservation ability on these special cases, as illustrated in Fig. 6.
We apply the scale normalization as proposed by Lindeberg [44].The final image reconstruction from the multiscale filtered orientation scores is obtained via .N o } and where N o and S represent resp.the number of orientations and the set of spatial scalings, and σ s ,σ o η,norm denotes the left-invariant Gaussian matching filters on the orientation score [65].The maximum filter response is calculated over all orientations per position.The method is validated on the public databases DRIVE and STARE and is computationally efficient.With a sensitivity and specificity of 0.7744 resp.0.9708 on DRIVE, and 0.7940 resp.0.9707 on STARE, the proposed algorithm outperforms most current segmentation schemes.Interestingly, it can also deal with complicated vessel configurations.For more details see [65,66].

Multi-scale and multi-orientation segmentation (BIMSO) for SLO images
For retinal images taken with a scanning laser ophthalmoscope (SLO) camera, new segmentation techniques are g-i the results of our method, and j-l the corresponding ground truth annotations by a human observer [57].From [65] required.Such cameras typically exploit two wavelengths, green and infrared and typically exhibit better contrast for hemoglobin (vessels and bleedings).However, they also exhibit speckle noise, due to which conventional analysis methods for RGB images often fail.Very few studies are dedicated to SLO image analysis (e.g., [63]).The proposed method [1], termed BIMSO (Brain-Inspired Multi-Scale and multi-Orientation), has four main steps: preprocessing, feature extraction, classification and postprocessing.The method is developed for the green channel (RGB and SLO).
Noise reduction is effectuated by a nonlinear gamma transform ( Ǔ f = α |U f | γ ).The orientation score is raised by a power factor γ > 1.The absolute value of the orientation score |U f | is taken because of the quadratic property of the cake wavelets, and α is determined by the sign Fig. 7 Top original SLO image (EasyScan, I-Optics Inc.).Bottom nonlinear orientation score gamma transformation with γ = 1.8.From [1] of the real part of the orientation score (Re(U f )).See Fig. 7.
After enhancement, a pattern recognition approach is taken to classify the vessels from the background.A normalized feature vector is constructed per pixel, including the intensity, the half-cake wavelet responses and 3 left-invariant multi-scale first order Gaussian derivatives and 8 secondorder Gaussian derivatives in the pre-processed orientation scores (two of the second-order derivatives are similar; therefore, there are only 8 unique second-order derivatives, so in total 11 Gaussian derivatives are used).The orientation score intrinsically allows for derivatives in the local directions (e ξ , e η , e θ ), and crossings are 'lifted' as different orientations in different orientation layers.The method is trained with a feed-forward neural network classifier and validated and tested on two different datasets (DRIVE for RGB and IOSTAR for SLO).
The BIMSO method detects blood vessels well, in particular the smaller ones in low contrast regions, and crossing vessels.The performance is better than the best supervised segmentation methods for RGB (e.g., [56]), for both the DRIVE and IOSTAR databases.For more details, see [1].

Crossing-preserving multi-scale vesselness
Frangi's vesselness filter [29], which is based on geometric relations between the principal curvatures (eigenvalues of the Hessian second-order matrix) to extract cylindrical shapes, has been extended to the orientation score domain, which makes it crossing-preserving [33].
In a single scale layer of the OS transform, the SE(2) vesselness is expressed as Frangi's vesselness formulation [29]: The coordinates for the calculation of the second-order derivatives can be defined with two possible orthogonal frames: the local moving frame of reference {∂ ξ , ∂ η , ∂ θ } determined by the tangent and normal of the orientation score geodesics, or the gauge frame {∂ a , ∂ b , ∂ c } [31] determined by the eigendirections of the Hessian matrix (H s,β U a f )(g) at scale s and g ∈ SE(2), normalized w.r.t. the β-metric.For the first we get: where the superscripts s,β indicate Gaussian derivatives at spatial scale s = 1 2 σ 2 s and angular scale 1 2 (βσ s ) 2 .The second frame leads to (with λ i the eigenvalues of the left-invariant Hessian): ), which can be interpreted as orientation confidence as defined by [31].The SE(2) multi-scale vesselness is computed by summation over all scales.Figure 8 shows that the multi-scale vesselness gives much better results at crossings and bifurcations.The gauge frame V a,b,c gives the best results.For more details, see [33].

Multi-orientation vessel tracking
Vessel tracking has the advantage over pixel classification that it guarantees connectedness of vessel segments.We have developed tracking via orientation scores [6], so classical difficulties as crossings, bifurcations, closely parallel vessels From [24] and vessels of varying width and high curvature can be dealt with naturally.
To handle bifurcations properly, one-sided kernels are exploited, constructed by decomposition of orientations scores in two opposite directions, weighted in the radial direction with an error function (cumulative Gaussian).
Two tracking methods are developed: Edge Tracking based on Orientation Scores (ETOS) and Centerline Tracking based on multi-scale Orientation Scores (CTOS).ETOS tracks both edges of the vessel simultaneously.The edge positions are detected in the orientation score, as a local minimum (left edge) and maximum (right edge) from the anti-symmetric imaginary part of the tangent plane.CTOS exploits the fact that the disturbing vessel central light reflex is filtered out by Gabor kernels with proper scales.The invertible orientation score has some distinct advantages over Gabor filtering: much lower computational costs as the orientation score filters for all scales simultaneously, whereas the Gabor filter only calculates a single spatial frequency at a time, and the orientation score filtering is more accurate.
Overall, ETOS outperforms CTOS.ETOS gives best results when applied on invertible orientation scores, can deal with many complex geometries, and gives reliable width measurements.

Sub-Riemannian geodesics in SE(2)
In the orientation score, vessels can also be tracked as geodesic curves [8,9,52].Geodesics (optimal curves minimizing some curve-length functional) are found by making proper use of the curved geometry of the domain.On the space of positions and orientations, a Riemannian geometry is defined, in which distances between tangent vectors are defined relative to the location in the score (e.g., via the left-invariant {∂ ξ , ∂ η , ∂ θ }-frame).Furthermore, since not all curves in the lifted domain SE(2) are natural some directions are prohibited, resulting in the definition of a sub-Riemannian geometry (compare this with the restricted motion of a car).Similarly, the paths of vessels are continuous and they do not exhibit sudden jumps to the left or the right.
The algorithms for sub-Riemannian geodesic extraction are implemented via an anisotropic fast marching scheme [49,52].As a result, the curve extraction procedure is both fast and robust.In summary, the sub-Riemannian geodesic extraction in SE(2) allows for the robust extraction of curves because: (1) crossing structures are disentangled in the orientation scores, and (2) high curvature (e.g., discrete jumps to the left or right, or sudden change of direction) is punished due to the restricted Sub-Riemannian geometry.See examples in Fig. 9.

Arterio-venous ratio (AVR)
The separation of vessels in arteries and veins is crucial, as their different physiology, flow and mechanical properties respond differently to disease development, e.g., arteriolar narrowing, a decrease of the artery calibers relatively to the vein calibers is an important early biomarker for diabetic retinopathy.Also tortuosity measures are expected to be different for arteries and veins [39].An auto-Fig.9 Data-adaptive sub-Riemannian geodesics give good and smooth tracking results, also over crossings.From [8] 123 mated artery-vein classification is required for high-volume screening.
The ratio of the arteriolar and venular diameters is called the arteriovenous ratio (AVR) and is classically computed from the six widest arteries and veins in a restricted zone around the optic disk [42].
We developed a novel method for artery/vein classification [25] based on local and contextual feature analysis of retinal vessels.Features are (a) the color, as arteries appear brighter and veins darker due to the oxygen content of the blood, (b) the transverse intensity profile, as arteries have a more pronounced central reflex to the camera flash, and (c) graph path properties of crossings and bifurcations of vessels, as these provide contextual information, because arteries never cross arteries and veins always cross arteries [17].
A non-submodular energy function is defined, integrating the features, and optimized exactly by means of graph cuts.The method was validated with a ground truth data set of 150 fundus images.An accuracy of 88.0 % was achieved for all vessels and 94.0 % for the 6 widest vessels.The contextual information especially benefits the classification of smaller vessels.For more details see [25].

Fractal dimension
The concept of fractal dimension was initially defined and developed in mathematics.It measures the complexity of self-similar objects that have the same patterns across different scales, e.g., trees and snowflakes.The vascular tree on the human retina also has self-similar branching patterns over different scales.Therefore, there is a growing interest in retinal fractal analysis, exploiting fractal dimension as a biomarker for discriminating healthy from diabetic retinopathy.
Fractal dimensions have been widely investigated.However, conflicting findings are found [3,12].This motivated us to especially investigate the stability and reproducibility of fractal dimension measurements.
The retinal vascular network is extracted by the method described in Sect.4.4.1.The region-of-interest (ROI) must be carefully specified for fractal dimension calculation, because it is a global measurement and different fundus cameras have different field-of-views.The ROI is determined via firstly locating the optic nerve head by the template-matchingbased method as described in Sect.4.3, and the optic disk radii are obtained.The fovea position is found at the global minimal intensity in a small region with 5*OD radii lateral distance to the optic nerve head.For the fovea resp.ONH-centered images, the circular ROI mask is centered at the fovea centralis, resp.ONH.For the fractal dimension calculation, we used the box counting method which paves the full retinal ROI image with square boxes (tiles) with different side-length (different scales) [48].In each box, various measurements are done: counting the number of boxes that overlap with the vessels, calculating the entropy of each box and calculating the probability of finding vessels near each pixel.Finally, these measurements are plotted against the box side-length on a log-log plot, and the fractal dimension is the slope of the fitting regression line.See Fig. 10.
We examined the stability of the fractal dimension measurements with respect to a range of variable factors: (1) different vessel annotations obtained from human observers, (2) different automatic segmentation methods, (3) different regions-of-interest, (4) different accuracy of vessel segmentation methods, and (5) different imaging modalities.Our results demonstrate that the relative errors for the measurement of fractal dimensions are significant and vary considerably according to the image quality, modality and the technique used for measuring it.So automated and semiautomated methods for the measurement of fractal dimension are not stable enough, which makes fractal dimension not a proper biomarker in quantitative clinical applications [35].

Bifurcations/crossings detection
We have developed a fully automatic bifurcation and crossing detection algorithm called Biologically Inspired CRossing detecting in Orientation Scores (BICROS) [55], see Fig. 11, which does not depend on vessel segmentations.The precision of the junction detections is high, and the method can discriminate with high accuracy between crossings and bifurcations.Interestingly, junctions of small vessels, which are typically missing in vessel segmentations, are well detected by BICROS (Fig. 12).Through this, the proposed hybrid method outperforms state-of-the-art techniques [5].It performs well on both RGB color and noisy SLO retinal images.It provides the orientation-augmented landmarks for image stitching, and because the orientations and the tracks of the branches are provided, useful other biomarkers can be extracted such as segment lengths and bifurcation angles.

Curvature/tortuosity
Vessel curvature is an important biomarker [13,34,39,53].Most methods require pre-segmentation and center line extraction, but our brain-inspired orientation score method [11] works directly on the image data by fitting so-called exponential curves in the orientation score.An exponential curve is a curve whose spatial projection has a constant curvature (Fig. 13).
The curvature values are computed directly from tangent vectors of exponential curves that locally best fit the data: The coefficients of the tangent vector c ξ , c η , c θ are expressed in the moving frame (see Eq. ( 1), [24,30]).Principal directions are calculated from the eigenvectors of the Hessian matrix, but now in a curved domain, for which we need left-invariant derivatives: As left-invariant derivatives do not commutate, e.g., ∂ θ ∂ ξ U = ∂ ξ ∂ θ U , we need to make the asymmetric Hessian matrix symmetric via H μ U = M μ (HU ) T M μ 2 (HU )M μ , so we can calculate the eigensystem of H μ U .The diagonal matrix M μ −2 compensates for the difference in units between the orientation and spatial dimensions, for details see [11].
A confidence measure is acquired by calculating the Gaussian Laplacian in the plane orthogonal to the tangent direction of the vascular geodesic.From all pixelwise curvature measures, several global measures are derived, as the mean and standard deviation.The method was validated on synthetic and retinal images.

Next steps
The computer-aided diagnosis algorithms are mostly written in Mathematica 10 (Wolfram Research Inc.).Processing times are 1-15 min per image, but this is shortened by using multi-core (32) servers, both in China as in the Netherlands, and developing GPU-based implementations, in particular for the OS-based algorithms.
A dedicated workstation has been designed, to program any pipeline of processes, in batch mode, and with automated report generation (see Fig. 15, see also [19]).Data are acquired with different fundus cameras (EasyScan, DRS, Topcon).The workstation is also used for annotations: at least three experienced ophthalmology experts diagnose and annotate the bio-markers on the fundus image datasets as ground truth or reference standard.We have established a detailed workflow protocol to minimize variations.
The crucial validation phase, not reported in this, paper, is ongoing.The algorithms are currently evaluated in the collaborating hospitals: image and extensive metadata have been acquired of 3000 diabetes patients in Shengjing Hospital, and over 20.000 normal images with first expert reading are acquired at He Vision shops in Shenyang, China.Studies are ongoing to correlate the rich set of geometric biomarkers described in this paper with the extensive Chinese metadata obtained in the clinical setting, and find the most predictive and effective (combination of) geometric biomarkers extracted from the retinal fundus images for DR, glaucoma, and AMD.

Conclusion
The brain-inspired multi-orientation approach turned out to be highly successful in dealing with many complex vessel geometries, such as crossings, bifurcations, curvature analysis without segmentation, enhancement and segmentation.The theory for orientation scores in SE( 2) is now well developed.The generalized approach of a Lie-group analysis by the visual front-end is currently a broad topic of research.The paper gives a review of many current geometric developments for retinal vessel analysis with the goal of establishing biomarkers for screening for early signs of a range of systemic, neurodegenerative and cardiovascular diseases.
In the parallel deep machine learning field major developments take place.The RetinaCheck team ended at 17th position in the recent Kaggle challenge on Diabetic Retinopathy classification (https://www.kaggle.com/c/diabetic-retinopathy-detection), of 661 participating teams.It can be foreseen that the merge of geometric analysis and reasoning with deep learning, based on the large volumes of data acquired in this project, can be highly successful.This is ongoing development, and will be reported in forthcoming papers.

Fig. 1
Fig. 1 Clinical data acquisition at the Department of Endocrinology of Shengjing Hospital, Shenyang, China.The camera is an EasyScan scanning laser fundus camera (i-Optics Inc.The Hague, the Netherlands).Extensive clinical metadata are co-recorded for validation

Fig. 3 a
Fig. 3 a Real and b imaginary part of the cake kernel in the spatial domain, c Fourier contours at 70 % of the maximum for all orientations and d B(ρ) with ρ = 0.8ρ n , the Nyquist frequency ρ n and N = 8.From[6]

2σ 2 and 2 s
where σ s and σ o are used to define the spatial scale 1 2 σ and orientation scale1  2 σ 2 o of the Gaussian kernel.Note that

Fig. 4
Fig.4 Exemplary retinal image patches used for template training.Top row positive patches.Bottom row negative patches.From[7]

Fig. 8 a
Fig. 8 a Retinal image f and multi-scale vesselness filtering results for b the Frangi filter V Fr ( f ) and our two methods c V ξη ( f ) resp.d V a,b,c ( f ) (left to right).From[24]

Fig. 10
Fig.10 The pipeline for calculating the fractal dimension from a color fundus image.From[35]

Fig. 11
Fig. 11 Top row pipeline of the BICROS method.a input image, b preprocessing result, c construction of the orientation score, d candidate selection, and e feature extraction and classification into bifurcations (orange) and crossings (blue).Bottom row The orientation column at