Euler Characteristic Surfaces
Gabriele Beltramo, Rayna Andreeva, Ylenia Giarratano, Miguel O. Bernabeu, Rik Sarkar, Primoz Skraba
EE ULER C HARACTERISTIC S URFACES
A P
REPRINT
Gabriele Beltramo ∗ Rayna Andreeva Ylenia Giarratano Miguel O. Bernabeu Rik Sarkar Primoz Skraba School of Mathematical Sciences, Queen Mary University of London, London, E1 4NS, UK [email protected],[email protected] School of Informatics, University of Edinburgh, Edinburgh, EH8 9AB, UK [email protected], [email protected] Usher Institute, University of Edinburgh, Edinburgh, EH16 4UX, UK [email protected], [email protected]
February 17, 2021 A BSTRACT
We study the use of the Euler characteristic for multiparameter topological data analysis. Eulercharacteristic is a classical, well-understood topological invariant that has appeared in numerousapplications, including in the context of random fields. The goal of this paper is to present theextension of using the Euler characteristic in higher-dimensional parameter spaces. While topologicaldata analysis of higher-dimensional parameter spaces using stronger invariants such as homologycontinues to be the subject of intense research, Euler characteristic is more manageable theoreticallyand computationally, and this analysis can be seen as an important intermediary step in multi-parameter topological data analysis. We show the usefulness of the techniques using artificiallygenerated examples, and a real-world application of detecting diabetic retinopathy in retinal images. K eywords Topological Data analysis · Filtrations and Bi-filtrations · Euler Characteristic Curves and Surfaces · Imageclassification
The field of topological data analysis (TDA) has attracted a lot of research over the last few years. I this field, the mostwidely used tool is persistent homology, which provides a stable summary of a space/dataset over an entire range ofparameter choices/scales. The application of this technique has been largely limited to one dimensional parameterspaces due to both theoretical and computational challenges.In this paper, we present an alternative approach to studying data where higher-dimensional parameter spaces arenaturally present. Rather than using homology, we opt for a simpler topological invariant – the Euler characteristic.We obtain highly efficient algorithms to compute summaries of datasets over multiple parameters. We refer to thesummaries as
Euler characteristic surfaces or Euler surfaces for short. The goal of this paper is to illustrate that Eulersurfaces can provide insight into the data over multidimensional parameter spaces .The Euler characteristic χ makes an appearance in many different fields of topology and geometry including: algebraictopology [1], differential geometry [2], and stochastic geometry and topology [3]. It has been generalized to highlyabstract settings such as enriched categories [4], and can be seen as generalized measure [5]. One of the trulyremarkable aspects of the Euler characteristic is that it has allows for a local description which directly enables its ∗ First author contacts. Email: [email protected]; Personal webpage: https://gbeltramo.github.io/ a r X i v : . [ m a t h . A T ] F e b PREPRINT - F
EBRUARY
17, 2021efficient computation. In this paper, we provide a general algorithm for two common settings in TDA and a Pythonpackage.The idea of topological invariant of a sequence of spaces – e.g. as found with persistent homology – appears in the Eulercharacteristic domain as the
Euler characteristic curve (ECC), and has been used for topological inference in a rangeof applicaitons [6]. We provide a brief overview of this work in the following section, with a focus on a data-drivenapproach, and efficacy in real datasets.Our main contribution is to investigate the multi-parameter setting. Beginning with the work by Carlsson andZomorodian [7], there has been a large number of approaches proposed to deal with multi-parameter persistence.Approaches include directions such as the rank invariant [7, 8, 9], microlocal analysis [10], and higher-dimensionalanalogues of persistence diagrams [11, 12]. These all capture related but somewhat different concepts. To the best ofour knowledge, the rank invariant is the only case where implementations exist, and although they perform well, thealgorithms do not scale in the same way as one-dimensional persistence. While multidimensional persistent homologyremains an active area of research, we note that there is substantial evidence that in many settings – particularlyin presence of randomness – there is a surprisingly little loss of information in going from homology to the Eulercharacteristic. It has been observed that in many random models, at any given parameter, the homology of a singledimension dominates [13, 14]. Therefore, studying the Euler characteristic at different parameter values can provide alot of topological information.As illustrated by the applications of the ECC, the Euler characteristic provides a useful functional summary of data,which can readily be used for classification – especially as closed form expressions often exist. The extension tohigher-dimensional parameter spaces clearly enables more discriminative summaries. More importantly, by consideringthe difference of Euler characteristics, we can identify interesting regions of parameter space .We present several instances where the “shape" of Euler surfaces and the difference of Euler surfaces provide interestinginformation about the underlying generating process and the parameter space. In addition to simulated data, we alsopresent a real-world application: detecting diabetic retinopathy (abbreviated DR). In diabetic patients, one of the earlymanifestations of this disease is change of the structure of blood vessels in the retina (See Figure 13). Accurate detectionof such changes may help in early detection of the disease and prevention of significant damage. Recent work inthis area has been a series of machine learning approaches aimed at automated detection of DR and other diseasesfrom retinal images [15]. These works have largely used methods such as neural networks, which are accurate onlywith large training data volumes and are not easily interpretable. In diagnostic medicine, datasets are often small andinterpretability is paramount. On multiple datasets of retinal images, we illustrate two key points. First: the ECCis already effective at detecting DR, and secondly: expanding to the multiple parameters can yield insights into thequalitative differences between the blood vessels in healthy patients and those suffering from DR.The overall goal of this paper is to highlight this natural extension to existing work, as we believe this may lead toan important scalable, multi-parameter technique to the TDA toolbox. The structure of the paper follows our maincontributions.• We define the Euler surface corresponding to bi-filtrations and higher-dimensional parameter spaces, and relateto current directions in TDA (Section 2);• Give efficient algorithms for a variety of input including cubical and simplicial complexes arising fromembedded point clouds along with the a Python package for computing the Euler surfaces (Section 5);• Show that Euler curves and surfaces are useful for a variety of classification tasks both for data generated byrandom models (Section 6), as well as real-world medical data (retinal images) (Section 7);• Most importantly show how Euler surfaces can give insight into the structure of datasets by highlighting“interesting" areas of the parameter space (Section 7).
This work can be seen as an attempt to understand multi-parameter filtrations, which avoids the difficulties inherent inmulti-parameter persistence by considering a simpler topological invariant: namely the Euler characteristic. We do notrecount the numerous approaches to multi-parameter persistence here as it is not directly applicable to this work andwill be obvious to experts in the field, while we certainly lose quite a bit of information in using this simpler invariant,we gain a readily computable and applicable approach to data analysis.The connection of Euler characteristics and data analysis goes back to the Kac-Rice formula which gives the expectednumber of critical points of a sufficiently nice random field [3]. This is most naturally thought of as the study of arandom function on a space. Taking the sublevel/superlevel sets of the random function yields a one-dimensional2
PREPRINT - F
EBRUARY
17, 2021filtration, which in turn for every function gives a piecewise constant integer valued function – this is called the Eulercharacteristic curve. As the input is random, it is natural to take the expectation. Due to what can be understated asfortuitous, there exists a closed-form formula for the expected Euler characteristic curve, which is called the Gaussiankinematic formula. This applies to wide range of spaces, see [3] for an in-depth account.This has been applied to fMRIs [16], cosmology [17], and more recently various machine learning classificationproblems [18]. There has also been research into efficient streaming algorithms for their computation on image data[19]. In terms of multi-parameter settings, we mention [20], which proves certain convergence properties of a Eulersurface arising from smoothing a Gaussian random field (GRF).Multiple Euler characteristic curves have also been used for shape analysis in the form of the Euler characteristictransform (ECT) [21, 22, 23]. In Section 4, we comment and describe the relationship between the two constructions.This can be seen indirectly as an alternative approach to multi-parameter settings, with both applications and interestingtheoretical properties. The ECT is closely related to the Euler integration [24, 25, 5], which like the Euler charateristichas appeared in several different mathematical areas.
We begin with the definition of the Euler characteristic:
Definition 2.1.
Let X be a CW complex of dimension k . The Euler characteristic χ ( X ) of X is χ ( X ) = k (cid:88) i =0 ( − i n i , (1) where n i = |{ σ ∈ X | σ is i dimensional }| , or the number of i -dimensional cells. There are many other ways to characterize the Euler characteristic, including the alternating sum of Betti numbers orthe integral of curvature over a Riemannian manifold (or appropriate triangulation). As our primary interest is (finite)data, we do not further recount the basic properties of the Euler characteristic, referring the reader to any standardtext on algebraic topology [1] or differential topology [2]. Furthermore, we may assume that we are generally in asetting where different notions of Euler characteristic coincide. While our definition is in terms of CW complexes, wefocus on two special cases which represent constructions which are either direct representations of data or are commonconstructions from data. That is, we restrict ourselves to1. simplicial complexes2. cubical complexes/imagesIn the case of simplicial complexes, we consider primarily proximity based complexes such as the alpha complex andthe Vietoris-Rips complex. While we give different constructions, we generally denote the underlying complex by K . Remark 1.
While we use the corresponding structure of these special cases in the algorithms, there is a clearmodification for more general complexes, e.g. cellular complexes, which may arise from other types of preprocessingsuch as collapses using Discrete Morse Theory. However, as the algorithms are linear in the size of simplicial complexes,it is likely that any additional preprocessing is likely to increase computation time, so we do not consider it here.
We begin with the one-dimensional parameter case. Due to its connection with persistent homology, it is the mostwell-studied and familiar case in topological data analysis (TDA). The basic object of study is no longer one space, buta sequence of spaces called a filtration . A filtration is a increasing sequence of nested spaces: ∅ ⊆ X ⊆ X ⊆ · · · ⊆ X m . (2)The indexing set may be discrete ( Z ) or continuous ( R ). One of the most common ways to construct filtrations in TDAis as sub-level sets of functions. Given a real-valued function f : X → R and a threshold α , we may obtain a space X α = f − (( −∞ , α ]) . By varying α , we obtain the sub-level set filtration (resp. the super-level set filtration induced by f − ([ α, ∞ ) ) inducedby f . In many cases of interest, such as piecewise-linear functions on simplicial complexes, these filtrations aretopologically equivalent to a filtration induced by a function which is piecewise constant on each simplex. Namely for acell, σ ∈ K , we have a function g ( σ ) = max x ∈ σ f ( x ) PREPRINT - F
EBRUARY
17, 2021 (a) (b) (c) (d)
Figure 1: From left to right: × gray-scale image, full cubical complex Q of the image in (a), finite point set in R ,and Delaunay complex D of the points in (c). (a) Filtration of the gray-scale image in Figure 1a.(b) Filtration of the finite set of points in Figure 1c. Figure 2: Filtrations of the example data in Figure 1. In both (a) and (b) the filtration parameters are displayed aboveeach subcomplex. The last subcomplexes in the sequence are the full cubical complex Q and the Delaunay complex D .Under this definition, the sublevel sets of the function forma filtration by subcomplexes since g ( τ ) ≤ g ( σ ) for every τ, σ ∈ K with τ a face of σ .Since our main interest lies in applications, we restrict ourselves to this piecewise constant setting, with an increasingsequence of finite complexes. In this case, we may assume without loss of generality, that the indexing set is discrete,although for exposition, we refer to the function value of a simplex as to when it enters the filtration.When these restriction is valid has been well-studied, but we note that this is the case well-behaved functions of finitecomplexes. For the interested reader, we point out that the classes of functions which are well-behaved include Morseand constructible functions [26, 27], but in many cases even these conditions may be relaxed.4 PREPRINT - F
EBRUARY
17, 2021 (a) (b)
Figure 3: The plot in (a) is the Euler characteristic curve of the image in Figure 1a, where sublevel sets are on valuesbetween and . The plot in (b) is the Euler characteristic curve of the points in Figure 1c, where the sublevel setsare on values a s such that there is a one simplex difference between Q a s − and Q a s for each ≤ s ≤ m .We make one further simplifying assumption: that cells enter the filtration one at a time. This is neither necessarynor restrictive but is rather done for the exposition of the algorithms. A simple condition which ensures this, is theassumption that all cells are assigned a unique function value.We refer to a piecewise constant function which induces a filtration as a filtering function . We now define the Eulercharacteristic curve. Definition 2.2.
Given a filtering function h : K → R , the Euler characteristic curve induced by h is the integer-valuedfunction E h : R → Z defined by E h ( α ) = χ ( f − (( −∞ , α ])) , (3) for each α ∈ R . Examples of the Euler characteristic curves of are given in Figure 1 are in Figure 3. As described in Section 1.1, thiscurve has been used in a number of difffernt applications.
Our goal is to extend this to multi-parameter filtrations. In generalizing, we consider filtrations which are Cartesianproducts of one parameter filtrations.
Definition 3.1. A k -parameter filtration is the Cartesian product of k one-parameter filtrations. Our restriction to products may seem restrictive. However, any finite poset may be embedded into a product of linearorders, which is referred to as the order dimension or the Dushnik–Miller dimension. While in general, computing theminimal embedding dimension of of a poset is NP-hard, there are many special cases which are known. Constructivetechniques for embedding posets into Cartesian products is interesting but we leave it for further work, as our maininterest in applications come from Cartesian products induced by functions in R d (although there are many cases ofinterest which do not fall into this category). In the interest of readability, we focus primarily the case of bi-filtrations(also since these are also readily visualised). Definition 3.2.
Let h : K → R , h : K → R be two filtering functions on a complex, and h = ( h , h ) : K → R the function defined by h ( σ ) = ( h ( σ ) , h ( σ )) for each σ ∈ K . Given two sets of monotonically increasing values R = { a s } m s =0 , R = { b t } m t =0 and defined K s,t = h − (( −∞ , a s ] × ( −∞ , b t ]) , the sublevel set bi-filtration of K induced by h on R , R is the grid of nested subcomplexes K , ⊆ K , ⊆ · · · ⊆ K ,m ⊆ ⊆ ⊆ K , ⊆ K , ⊆ · · · ⊆ K ,m ⊆ ⊆ ⊆ ... ... . . . ... ⊆ ⊆ ⊆ K m , ⊆ K m , ⊆ · · · ⊆ K m ,m (4) We say that h : K → R is a bi-filtering function on K . PREPRINT - F
EBRUARY
17, 2021
Definition 3.3.
Let h : K → R be a bi-filtering function on K and R = { a s } m s =0 , R = { b t } m t =0 two set ofmonotonically increasing real values. The Euler characteristic surface induced by h on R , R is the matrix of integervalues S h = χ ( K , ) , χ ( K , ) , · · · χ ( K ,m ) χ ( K , ) , χ ( K , ) , · · · χ ( K ,m ) ... ... . . . ... χ ( K m , ) , χ ( K m , ) , · · · χ ( K m ,m ) , (5) where K s,t = h − (( −∞ , a s ] × ( −∞ , b t ]) . The last column and last row of S h are the Euler characteristic curves of h and h respectively. So the Eulercharacteristic surface contains all the topological information of E h and E h at the values in R and R , plus theinformation coming from intersection of sublevel sets of h and h .In practice, we have two forms of data we consider: digital images, which are treated as a form of cubical complex, andsimplicial complexes, where we focus on proximity complexes on point clouds. Definition 3.4.
Let M be a n -by- n two-dimensional gray-scale image. The cubical complex Q of M is defined bythe union of squares C i,j = [ i, i + 1] × [ j, j + 1] and their subfaces for each pixel ( i, j ) of M . We define the pixelintensity filtering function h M : Q → R setting(i) h M ( C i,j ) = v i,j for each ≤ i ≤ n and ≤ j ≤ n , where v i,j is the intensity of pixel ( i, j ) ;(ii) h M ( C ) = min C ⊆ C i,j h M ( C i,j ) for each element C of Q .In the same way, given a three-dimensional gray-scale image M , h M is defined by voxels intensities v i,j,k setting h M ( C i,j,k ) = v i,j,k and h M ( C ) = min C ⊆ C i,j,k h M ( C i,j,k ) for each voxel ( i, j, k ) and each element C in Q . We also consider simplicial complexes built on point sets in R d . We focus on proximity complexes – complexes whichdepend on distances between the points. These include common constructions in TDA such as ˇCech complexes,Vietoris-Rips complexes, alpha and Delaunay complexes [28].Again we consider the image M in Figure 1a, and the finite set of points X in Figure 1c. The filtrations induced on thisdata by the filtering functions h M and h X are in Figure 2a and Figure 2b respectively. These functions may be based ondistance or some other function(s). In our experiments, we use the Delauanay complex as most of our experiments arein two or three-dimensional space. Remark 2.
A key advantage to using the Euler characteristic is that it is defined pointwise. That is, the value isdetermined by each complex, the maps between the spaces have no effect on the invariant, i.e. different maps will leadto the same Euler characteristic surface. While this does make it a weaker invariant, it also implies that the problemspresent in multidimensional persistence are avoided.
Example: Euler characteristic surfaces of random images.
We show with an example that the Euler characteristicsurfaces of pairs of images can contain useful information for distinguishing between different classes in a dataset,while the Euler characteristic curves of the same images do not.Our data consists of families of pairs of n × n random images M , M , whose pixel intensities are uniformlydistributed and have an expected correlation p . These are generated by drawing three sample values for each pixel ( i, j ) .In particular, we draw x, v , v from independent uniform distributions U (0 , , U (0 , , U (0 , . If x ≤ p , weset M [ i ][ j ] = M [ i ][ j ] = (cid:98) v (cid:99) . Otherwise, we set M [ i ][ j ] = (cid:98) v (cid:99) and M [ i ][ j ] = (cid:98) v (cid:99) . Hence, we set M [ i ][ j ] and M [ i ][ j ] to the same random integer with probability p and to independently drawn random integers with probability (1 − p ) .Given a pair of gray-scale images M , M , we obtain an Euler characteristic surface by bi-filtering on the pixelintensity filtering functions h M , h M . In this setting, we can derive the expected value of the Euler characteristic χ ( Q s,t ) = χ (cid:0) h − (( −∞ , s ] × ( −∞ , t ]) (cid:1) for each pair of thresholds ≤ s ≤ m , ≤ t ≤ m . We know that a vertex ( i, j ) is in Q s,t if and only if at least one of the squares that include it is in Q s,t . The same holds for edges in Q s,t .Thus, given the probability of having squares in Q s,t , the probabilities of having vertices and edges can be derived. By6 PREPRINT - F
EBRUARY
17, 2021 (a) (b) (c)
Figure 4: Contour plot of S ( h M ,h M ) of a pair of random images M , M with correlation coefficient equal to . in(a) and equal to . in (b). The difference between these two Euler surfaces is the Euler terrain in (c).the definition of M and M , in terms of random values sampled from uniform distributions, it follows that P ( C i,j ∈ Q s,t ) = P (cid:16) h M ( C i,j ) < s and h M ( C i,j ) < t , with h M ( C i,j ) = h M ( C i,j ) (cid:17) · p + P (cid:16) h M ( C i,j ) < s and h M ( C i,j ) < t , with h M ( C i,j ) , h M ( C i,j ) independent (cid:17) · (1 − p )= min { s, t } · p + s · t · (1 − p ) (6)where ≤ s, t ≤ and C i,j is any square in the two-dimensional cubical complex Q . Then, because thevalues of different pixels are independent of each other, the probability that a vertex/edge σ (cid:48) belongs to Q s,t is − (1 − P ( C i,j ∈ Q s,t ) k ) , where k is the number of squares containing σ (cid:48) .Finally, there are n · n squares in the cubical complex Q of M and M . These contain ( n + 1) · ( n + 1) vertices, subdivided into ( n − · ( n − internal vertices contained into squares each, n −
1) + 2( n − boundary vertices contained into squares each, and corner vertices contained in a single square. Moreover, there are n ( n + 1) + n ( n + 1) − n − n interval edges contained in squares each, and n + 2 n boundary edgescontained in a single square. Hence the expected value of χ ( Q s,t ) is E [ χ ( Q s,t )] = ( n − · ( n − · [1 − (1 − P ( C i,j ∈ Q s,t ) )]+ ( n · ( n + 1) + n · ( n + 1) − · [1 − (1 − P ( C i,j ∈ Q s,t ) )]+ ( n · n + 2 n + 2 n + 4) · P ( C i,j ∈ Q s,t ) . (7)Figure 4 displays the contour plots of expected Euler characteristic surfaces for two different values of p . The imagesdefining them have the same uniformly random distribution of pixel values, so the corresponding expected Eulercharacteristic curves are all equal. However, the Euler characteristic surfaces are able to discriminate between pairsof images with different correlation p , see the Euler terrain in Figure 4c. Furthermore, it is readily apparent thatthe structure of the of surfaces is encoding information about the correlation. Since this is a toy model, we do notcharacterize the behaviour rigorously here. However, it does serve as an indication as to the structural informationwhich is encoded in Euler surfaces. This can be particularly useful for comparing models with real-world data, as theshape of the surface may indicate what behavior in data, a model is failing to capture. This is especially intriguing asclosed form expressions often exist Euler characteristics [3, 14] or alternatively, as a well-behaved functional, it can beempirically estimated. Open Problem.
How Euler surfaces for models with multiple parameters differ for varying values of parameters?
Difference of Euler surfaces.
Expected Euler surfaces are most often well-behaved with closed-form solutions oftenexisting. Alternatively the empirical expected Euler characteristic can be efficiently computed as the surfaces are definedpointwise. Therefore, when comparing Euler surfaces which come from two classes, a natural object to investigateis the difference of average Euler surfaces . While this is computed pointwise, the resulting structure of this objectprovides insight into parameter values which differentiate between the two classes.
Definition 3.5.
Given two sets of Euler surfaces over a common parameter space denoted { K s,t } and { K (cid:48) s,t } , theirdifference at point s, t is defined as ECS(
K, K (cid:48) )[ s ][ t ] := 1 |{ K s,t }| (cid:88) χ ( K s,t ) − |{ K (cid:48) s,t }| (cid:88) χ ( K (cid:48) s,t ) PREPRINT - F
EBRUARY
17, 2021For brevity, we refer to this object as the
Euler terrain or simply as the terrain . One important point for this objectis that if we are comparing two distinct random processes where the expected Euler characteristic exists and thecorresponding surface is continuous (after scaling), replacing the empirical averages with expectations results in theexpected difference. By linearity of expectations, the terrain also exists, is continuous, and provides a comparison oftwo processes.In experimental instances where the number of samples for the two classes is highly unbalanced (or the number ofsamples is small), it makes sense to further normalize by the standard deviation of each Euler characteristic. We refer tothis as a normalized terrain , given by
ECS(
K, K (cid:48) )[ s ][ t ] := ECS ( K, K (cid:48) )[ s ][ t ]sd( χ ( K )[ s ][ t ]) + sd( χ ( K (cid:48) )[ s ][ t ]) where sd is the standard deviation. Terrains provide interesting information about the data which we investigateexperimentally in Sections 6 and 7. The Euler Characteristic Transform (ECT) [21, 22, 23] was originally considered for three-dimensional shape analysis.Here we present the basic form and refer the reader to [22] for more general and in-depth presentation.Given a compact shape in R d , consider the space of directions given as points on S d − . For each direction consider thesub-level set filtration given by the height function in that direction. This yields an Euler Characteristic Curve for eachdirection. There is the surprising theorem: Theorem 4.1 (Theorem 3.4 [22]) . If M and M (cid:48) are two constructible subsets of R d , then equivalence of the Eulercharacteristic transforms implies that the sets are equal, i.e. ECT( M ) = ECT( M (cid:48) ) ⇒ M = M (cid:48) This theorem essentially implies that we do not lose any information when passing to a collection of Euler characteristiccurves. There has been substantial interest toward understanding how finite sampling affects reconstruction.A natural question is how do Euler surfaces relate to the ECT. First, the ECT uses many filtering functions. In itstheoretical set-up, for a shape in R d , a function is constructed for each point in S d − . Even in the approximate case,many different directions are needed. As we are primarily interested in visualizing the results, we restrict to only a fewfunctions (most often only two). Hence the first (rather obvious) connection is that from an ECT, the choice of any k directions will result in a k -dimensional Euler Surface (or more accurately an Euler k -dimensional hypervolume).In our applications, there is often a natural real-valued function, i.e. pixel intensity, along with other parameters. Theabove theorem also applies to this setting. Corollary 1.
Given a function on a subset K ⊂ R d , the graph of the function can be reconstructed by the ECT.Proof. Since the domain of the function is a subset of R d , the graph can be embedded as a shape in R d +1 . The abovetheorem then immediately implies the result.One observation we make is that in the above case, we need not take all directions on S d . Rather, we can consider eachlevelset of the function separately. Each levelset is a set in R d , which can be reconstructed by considering directionsin S d − . Note that this does not reduce the number of directions considered but can be more intuitive when we arestudying real-valued functions.A second observation is that this argument can be iterated, allowing us to reduce to the case of reconstructing shapes in R with directions in S . While we believe this reduction can be useful in further analysis of the ECT, in this setting,the Euler Surfaces we use are a point in S for each function (see Figure 5). This perhaps best illustrates the connection,where the Euler Surfaces are a sampling of the ECT. As we will show, despite not characterising the shape completelyas the ECT, the Euler surfaces can still show us useful information, so we present the following open problem: Open Problem.
Given an ECT, is possible to automatically choose some set of “interesting" Euler surfaces which canhelp us understand the underlying structure?
In addition to being highly useful for data analysis, this can provide some better insight into the theory of ECTs.8
PREPRINT - F
EBRUARY
17, 2021Figure 5: The realtionship between the ECT and the ECS for functions on R . For every direction α , the ECT computesthe Euler characteristic curve based on the filtration arising from the height function in that direction. In the standardsetup, it would be a direction on the sphere S , but it could equivalently be all “directions” in S × R . The 2-dimensionalEuler surface we consider is a direction and the function sublevel sets, so the Euler surface can be thought of asrepresenting the slice shown. In this section, we describe efficient algorithms for the computation of Euler characteristic surfaces of image and pointdata. An implementation of these is provided by the euchar Python package, the source code of which is available onGithub at https://github.com/gbeltramo/euchar.
First, we describe an algorithm for the computation of the Euler characteristic surface of a pair of gray-scale images M , M . In particular, Algorithm 1 returns the surface of the bi-filtration of the sublevel sets of h : Q → R defined by h ( C ) = ( h M ( C ) , h M ( C )) for each C ∈ Q , where h M and h M are the pixel intensity filtering functions of M and M respectively. Discussion.
By definition of Cartesian product, we have that Q s,t = h − (cid:0) ( −∞ , a s ] × ( −∞ , b t ] (cid:1) is equivalent to Q s,t = h − M (cid:0) ( −∞ , a s ] (cid:1) ∩ h − M (cid:0) ( −∞ , b t ] (cid:1) . So each column of the Euler characteristic surface S h equals the Eulercharacteristic curve of h M with Q restricted to its top-dimensional cubes C such that h M ( C ) ≤ b t , because of theintersection with the cubical complex h − M (cid:0) ( −∞ , b t ] (cid:1) . Thus, a näive approach for computing the S h is to computeEuler characteristic curves multiple times with algorithms looping on the cubes in Q [19]. To improve over this,Algorithm 1 makes use of the following two strategies: (i) Precompute the possible Euler characteristic changes produced by adding a top-dimensional C into any Q s,t ,and use these to increase or decrease the values of S h ; (ii) Loop on each top-dimensional C only once, by modifying all columns of S h where C produces the samechange at the same time.In the following discussion, points (i) and (ii) above are shown to preserve the correctness of the naïve approachcomputing columns of S h independently.Using Euler characteristic changes as suggested in (i) is possible because the process of going from the empty abstractcubical complex to Q = Q , can be decomposed into steps at which a single ¯ C and its subfaces are added.9 PREPRINT - F
EBRUARY
17, 2021
Algorithm 2
Euler characteristic surface of bi-filtration on a pair of images.
Input: gray-scale images M , M , h : Q → [0 , m ] × [0 , m ] ⊆ R , and the pre-computed vector preCompChanges . Add a one pixel (voxel) thick outer layer to images, so that the new boundary pixels (voxels) are mappedby h into ( m + 1 , m + 1) S h ← ( m + 1) × ( m + 1) zeros matrix for each top-dimensional cube C in Q M do a s , b t ← h M ( C ) , h M ( C ) neigh , neigh ← h M , h M values in neighbourhood of C thresholds ← sorted values in neigh greater than b t , union m + 1 N C ← boolean matrix defined by ( neigh ≤ a s ) before C and ( neigh < a s ) after C for k = 1 to | thresholds | do N C ← boolean matrix defined by ( neigh ≤ thresholds [ k − N C ← element-wise AND of N C and N C l ← decimal integer of binary representation of N C for ˆ t = index of thresholds [ k − to index of thresholds [ k ] − do S h [ s ][ˆ t ] += preCompChanges [ l ] end for end for end for S h ← cumulative sum on columns of S h return S h This follows from the definition of the filtering functions h M and h M in terms of pixel (voxel) intensity values.Furthermore, at each such step, the change ∆ χ ¯ C in Euler characteristic of the current cubical complex is completelydetermined by the structure of elements adjacent to ¯ C . More precisely, defined the neighbourhood N ¯ C of ¯ C to be theset of cubes that intersect it, by Definition 2.1 ∆ χ ¯ C only depends on the numbers of cubes added into N ¯ C when ¯ C isadded.All possible Euler characteristic changes can be precomputed because there is a finite number of neighbourhoods N ¯ C . In particular, there are (3 d − such neighbourhoods in dimension d , meaning that there are Euler characteristicchanges to precompute for two-dimensional images and , , changes for three-dimensional images. For d = 4 ,the number of possible neighbourhoods is already a digits integer, making the computation and storage of theircorresponding changes impractical. Hence Equation (1) can be used to compute all the Euler characteristic changesfor d = 2 and d = 3 , which can then be stored in a vector preCompChanges using the binary representation ofneighbourhoods to index them. For example, consider the neighbourhood in Figure 6a corresponding to the binarymatrix (cid:32) (cid:33) , (8)and in turn to the binary sequence . Its Euler characteristic change is − and the decimal representation of itsbinary sequence . Thus − is stored as the -th element of preCompChanges .Point (ii) above is realized by the inner loop on lines − of Algorithm 1, where a s = h M ( ¯ C ) and b t = h M ( ¯ C ) sothat Q s,t is the first complex including ¯ C . The idea is to use preCompChanges to update the s -th row of S h at eachiteration. This can be done because Q s,t = h − M (cid:0) ( −∞ , a s ] (cid:1) ∩ h − M (cid:0) ( −∞ , b t ] (cid:1) , so χ ( Q s,t ) and χ ( Q s,t +1 ) can differ bya change ∆ χ ¯ C induced by ¯ C if and only if N ¯ C in Q s,t +1 has changed, i.e. if there is a top-dimensional cube C (cid:48) ∈ N ¯ C such that h M ( C (cid:48) ) = b t +1 . But all such changes depend on the h M values of top-dimensional cubes in N ¯ C greaterthan b t . Sorting and storing these in thresholds with m + 1 appended, it follows that the ranges of ˆ t -th columns of S h such that ˆ t is between two consecutive values of thresholds are such that the Euler characteristic change inducedby adding ¯ C is constant because N ¯ C does not change. So the elements of vector preCompChanges can be used online to update all ˆ t columns such that ˆ t ≥ j . For two-dimensional images, N ¯ C is a set of squares and their subfaces, while for three-dimensional images it is a set of cubes and their subfaces. PREPRINT - F
EBRUARY
17, 2021In conclusion, at the end of the loop on lines − , each entry S h [ s ][ t ] equals the change χ ( Q s,t ) − χ ( Q s − ,t ) , becauseall changes ∆ χ ¯ C induced by the top-dimensional ¯ C in Q s,t \ Q s − ,t have been considered. After the cumulative sumon columns of S h , it follows that S h [ s ][ t ] = (cid:16) χ ( Q ,t ) − χ ( ∅ ) (cid:17) + . . . + (cid:16) χ ( Q s,t ) − χ ( Q s − ,t ) (cid:17) = χ ( Q s,t ) − χ ( ∅ ) = χ ( Q s,t ) , (9)which is the required Euler characteristic surface entry.Ignoring the time required to precompute the vector of Euler characteristic changes preCompChanges , Algorithm1 has a worst case running time of O ( nm + m m ) , where n is the number of pixels (voxels) in M and M . Thisfollows because the inner loop on lines − takes O ( m ) operations in the worst case to update an entire row.However, compared to computing m + 1 Euler characteristic curves as proposed by the naïve approach at the beginningof this discussion, N ¯ C is computed only once for ranges of columns where it does not change, and entries S h areincremented and decremented without having to count subfaces of top-dimensional cubes in N ¯ C . (a) (b) Figure 6: Euler characteristic changes produced by adding a cube of maximal dimension in a two-dimensional cubicalcomplex Q s − ,t . In (a) the change is equal to − , while in (b) it is +1 . Given a finite point set X , which we assume being in general position, we provide Algorithm 3 for the computation ofthe Euler characteristic surfaces of a bi-filtration of a simplicial complex K built onto X . Moreover, it is assumed thatthis bi-filtration consists of sublevel sets of a h = ( h , h ) : K → R on monotonically increasing sets of real values R and R . Discussion.
In this case, when a simplex σ is added into a K s,t = h − (cid:0) ( −∞ , a s ] (cid:1) ∩ h − (cid:0) ( −∞ , b t ] (cid:1) its neighbourhooddoes not have a fixed structure. Thus it is not possible to precompute Euler characteristic changes as in Algorithm 1.However, if σ ∈ K s,t , then σ ∈ K s, ˆ t for each ˆ t ≥ t . So the change in Euler characteristic ( − dim ( σ ) , produced by Algorithm 3
Euler characteristic surface of bi-filtration on finite point set.
Input: abstract simplicial complex K , h = ( h , h ) : K → R , and sorted values in R and R . S h ← ( m + 1) × ( m + 1) zeros matrix for each simplex σ in K do v , v ← h ( σ ) , h ( σ ) a s , b t ← minimum values greater than v , v in R , R with binary search for ˆ j = j to m do S h [ s ][ˆ t ] ← ( − dim ( σ ) end for end for S h ← cumulative sum on columns of S h return S h PREPRINT - F
EBRUARY
17, 2021Table 1: Classification results obtained with Euler characteristic based feature vectors and logistic regression.
Outex_TC_00000
Features Avg test accuracy
Euler curve - pixel intensity . ± . %Euler surface - pixel intensity, Laplacian . ± . %MNIST Features Test accuracy
Euler curve - pixel intensity . %Euler surface - pixel intensity, top/bottom gradient . %adding σ into K s,t , also applies to K s, ˆ t for each ˆ t ≥ t . This property is used on line of Algorithm 3 to update the s -throw of S h for each σ . It follows that at the end of the loop on lines − each entry S h [ s ][ t ] equals χ ( K s,t ) − χ ( K s − ,t ) ,and the cumulative sum on columns of on line returns the desired Euler characteristic surface.Differently from Algorithm 1 (where value a s is mapped to index s , and b t to t ), it is necessary to find the indexes s, t of h ( σ ) and h ( σ ) within the sorted values of R and R . With a binary search this operation takes log ( m ) and log ( m ) operations for s and t respectively. Hence the worst-case running time of Algorithm 3 is O ( n ( log ( m ) + log ( m ) + m ) + m m ) , where n is the number of simplices in K .In the point set case, there is no inherent locality that can be exploited. In some cases, one could construct the complexlocally, but ultimately this is an improved technique of constructing the complex rather than any improvement incomputing the Euler characteristic, which counts the simplices as a sorted list. We conclude this section with a remark. Remark 3.
The algorithms presented here are quite straightforward and we include them primarily for completeness.This simplicity also leads to them being exceptionally efficient. One important open question is whether algorithms canbe made sublinear if we allow for approximations. This is particularly important for higher-dimensional parameterspaces, as the complexity rises exponentially in the dimension of the parameter space.
In this section, we first present the use of Euler surfaces for the classification for images and differenting betweenvarious random processes including random images and point processes. In Section 7, we look at a real-world dataset,but here we are able to investigate the resulting surfaces when the underlying process is known.
Here we consider two standard benchmarking image classification databases: the
OUTEX_TC_00000 test suite [29] andthe MNIST database of handwritten digits [30]. (a) (b)
Figure 7: The images in (a) are eight of the twenty-four patterns in the
OUTEX_TC_00000 test suite. The images in (b)are a random selection of the , handwritten digits from the MNIST database.We test the performance of Euler characteristic curves versus Euler characteristic surfaces as feature vectors of logisticregression. For this we use the implementation provided by the scikit-learn Python machine learning package [31],12
PREPRINT - F
EBRUARY
17, 2021 (a) (b)
Figure 8: Points sampled from two Clayton copula distributions. In (a) θ = 1 , while in (b) θ = 5 .together with a lbfgs solver. Each Euler characteristic curve is calculated directly on each gray-scale image M . Onthe other hand, to compute Euler characteristic surfaces we define a second image M for each one in Outex_TC_00000 and MNIST respectively. For the former dataset, we set M equal to the Laplacian of M , defined as the discreteconvolution of M with the kernel (cid:32) − − − − (cid:33) . Instead, for the MNIST dataset of images, we define M as the constant gray-scale image given by the top-down gradientof the same size of any MNIST image. Since MNIST images are × in size, we have that M [ i ][ j ] = (cid:98) · i (cid:99) .To obtain feature vectors for logistic regression, we subsample both Euler characteristic curves and surfaces bypreserving only one element in for curves, and one row and column in for surfaces. Also, we concatenate theresulting rows of Euler characteristic surfaces. Finally the components of the resulting feature vectors are normalizedto have mean and standard deviation . The results are given in Table 1, where we use average test accuracy as thescoring metric. In both cases feature vectors derived from Euler characteristic surfaces result in higher classificationscores. By combining multiple sources of information Euler surfaces outperform Euler curves. So, in setting inwhich data can be parameterized by more that one real-valued function, Euler surfaces are useful for encoding moreinformation than Euler curves into feature vectors employed in classification tasks. The aim of this experiment is to provide empirical evidence of the conclusions of the example given in Section 2. Therewe saw how the expected Euler characteristic surfaces of pairs of uniformly random images are able to distinguishbetween different levels of correlation p . Here, we synthetically generate correlated pairs of uniformly random three-dimensional images. To have a well defined notion of correlation we employ copula distributions [32]. In particular, wesample random points from a two-dimensional Clayton copula with uniform marginals U (0 , as implemented bythe copula R package, see [33]. An example is given in Figure 8, where , points were sampled from Claytoncopula distributions with different θ parameters.We make a pair of 3D random images M , M of shape n × n × n , by mapping the coordinate values of n pointssampled from a Clayton copula in R to the voxel intensities of M , M . Each θ results into a family of pairs of random3D images.We generate random pairs of 3D images with this method for both a Clayton copula distribution with θ = 1 andwith θ = 5 . The resulting Euler characteristic surfaces are averaged to obtain the surfaces in Figure 9a and Figure9b. The absolute difference of the two is in Figure 9c. As in the example in Section 2 the Euler characteristic curvesof any M and M are known to be have the same expected values. On the other hand, the Euler terrain in Figure 9cclearly shows that the average Euler surfaces of different families of images distinguish between different values of θ .13 PREPRINT - F
EBRUARY
17, 2021 (a) (b) (c)
Figure 9: Average Euler characteristic surfaces of random images generated with random points as in Figure 8. In (a)surface of images obtained with θ = 1 , and in (b) with θ = 5 . The absolute value of the difference of the surfaces in (a)and (b) is in (c). (a) (b) Figure 10: Point processes, corresponding to a Poisson point process in (a) and Hawkes cluster process in (b).
Here we study the properties of Euler characteristic surfaces of finite point sets. We start by providing an exampleusing Euler characteristic surfaces to encode the statistic properties of different point processes in R . In particular, weconsider (i) A Poisson point process in the unit square [0 , × [0 , of intensity λ = 400 ; (ii) A Hawkes cluster process, as described in Section of [34], whose cluster centers are generated as a Poissonprocess of intensity λ = 280 . We take a offspring intensity function ρ ( x, y ) = α πσ · e − σ ( x + y ) , with α = 0 . and σ = 0 . . Note that by the definition of Hawkes process and the fact that (1 − α )400 = 280 ,we have that the expected number of point generated by this process is .We generate finite point sets from both types of processes and define bi-filtering functions h = ( h , h ) : D → R on the Delaunay triangulation D of each of them setting• h = h X , the α -filtration values (which is related to the distance to D ) [35];14 PREPRINT - F
EBRUARY
17, 2021 (a) (b)(c) (d)
Figure 11: Average Euler characteristic surfaces obtained from the point processes in Figure 10. In (a) and (b) theaverage surfaces of points obtained from a Poisson and Hawkes cluster process respectively. In (c) the Euler terrainrepresenting their difference. In (d) the black area represent regions of the parameter space where the two averagesurfaces in (a) and (b) are significantly different.• h ( v ) = (cid:113)(cid:80) u ∈ U (cid:107) v − u (cid:107) | U | for each vertex v ∈ D , where U is the set of k -nearest neighbours of v , and h ( σ ) = max v ∈ σ h ( v ) for each σ ∈ D . This filtering function estimates the inverse of the density at eachvertex v , and extends it with the maximum to higher-dimensional simplices.The resulting average Euler characteristic surfaces are in Figure 11. Their absolute value difference in Figure 11c. InFigure 12, we show the complexes at built on two instances of the cluster process and two instances of the Poissonprocess, at parameters where the Euler surfaces differ and the difference can clearly be seen. We later use this sameapproach to understand a real data set. Here we present an application of the techniques to a real-world medical dataset on diabetic retinopathy. First, weintroduce the dataset and provide evidence that the Euler characteristic curve is useful for detecting the condition usingEuler characteristic curves (ECC). As the number of samples is small, the extension to surfaces is not the primarygoal. Rather, having established that the Euler characteristic does capture the relevant information, we come to oneof the main contributions of this paper. We use the surfaces to identify regions of parameter space which distinguishnormal samples from non-normal samples. This kind of visualization is highly useful, as it is simple to then look atthe underlying images and can be used to demonstrate the relevant direction/super or sub level set where differencebetween healthy and ill cases. 15
PREPRINT - F
EBRUARY
17, 2021 (a) (b) (c) (d)
Figure 12: Examples of the complex built at the radius and density corresponding to Region A in Figure 11(c) for acluster process (a,b) and the Poisson process (c,d)
Diabetic retinopathy (DR) is a common consequence of diabetes [36, 37], and a leading cause of blindness world-wide [38]. The disease causes a degeneration of blood vessels in the retina – Figure 13(a) and (c) shows a comparison.The signs of DR are intuitive to trained physicians, but not easy to detect computationally. We found that Eulercharacteristic based features are effective in detecting diabetic retinopathy on 2 different datasets (NHS Lothian andOCTAGON [39]) of size N = 51 and N = 43 respectively. As input, we have images of the blood vessels in the eyefrom different patients via a method called Optical Coherence Tomography Angiography. The images can be thought ofas a grayscale images with an intensity value where higher values correspond to the presence of blood vessels. In manyapplications, these blood vessels are first extracted as a graph or tree and then certain graph features are used as inputfor classification algorithms [40, 41]. OCTA.
Optical Coherence Tomography Angiography (OCTA) has been one of the most important prospective imagingmodalities in the retinal imaging domain over the last few years [42, 43]. It enables the visualisation of retinal bloodvessels in a rapid and non-invasive way [43]. Even though OCTA is not the established modality for diagnosing DRin the clinical practice, it has the potential to help doctors diagnose patients at the very early stages of DR, which arecrucial for prevention and treatment. There are early signs of DR which are more likely to be observed on OCTAimages before they are apparent on fundus examination [44]. Examples of OCTA images of healthy patients (Controls),patients with DR and patients with Diabetes who have not developed DR (NoDR) are provided in Figure 13.
Overview of experimental results with ECC.
We first compare the ECC between healthy and diseased patients on 2different OCTA datasets, coming from 2 different imaging devices (NHS Lothian (Optovue RTVue XR Avanti system(Optovue Inc., Fremont, CA)), N = 51 and OCTAGON (DRI OCT Triton system (Topcon Corp, Tokyo, Japan)), N = 43 ). The input images have resolution × and × for NHS Lothian and OCTAGON respectively.Using the ECC based on the intensity levels (sub and superlevel sets), we apply the method to 2 image classificationtasks: 2 class classification (Control vs. DR) and 3 class classification (Control vs. NoDR vs. DR). As seen in Table2 and 3, the ECC performs better than a baseline of 2 biomarkers and comparable to state-of-the-art approaches. Inparticular, we achieve AUC of . in the Control vs. DR study for NHS Lothian, . for OCTAGON and AUCof . for Control, . for NoDR and . for DR in the 3-class study (NHS Lothian). Most notably, we achieveaccuracy of for NoDR, which is the hardest class to classify due to the lack of well-defined signs in the image andthis accuracy is better than the transfer learning approach. Overview of experimental results with ECS.
Due to the small sample sizes, we do not attempt classification withEuler surfaces (as the curves already perform well). Rather, we use 2-parameter ECS to look for new topologicalfeatures. Using the same OCTA datasets as above, we use the image intensity as our first function. For the secondfunction, we used:1. the complement function as illustrated in 14, e.g. each point represents intensities in a range,2. the radial gradient image as seen in Figure 15b.The rationale for the first one is to identify useful intervals which identify the illness, where a range of intensitiesapproximately identify the blood vessel thicknesses, as thicker blood vessels correspond to higher intensities. Thesecond function is to identify if there is a distance from the center of the retina image (which mostly coincides with thecenter of the FAZ) where blood vessels begin to behave differently. These choices are not exhaustive but demonstrate16
PREPRINT - F
EBRUARY
17, 2021one of the main advantages to our approach – they are readily interpreted in the context of the data, i.e. they have aclear physiological interpretation.To test the approach, we construct two separate pipelines ( level set pipeline and radial gradient pipeline ), dependingon the function used. As seen in Figures 16, we identify regions of interest with resulting topological biomarkers foreach pipeline, for which we calculate a correlation coefficient with 2 known biomarkers. We detect that one of thesebiomarkers, vessel density, is strongly correlated with Pearson correlation coefficient of . with EC from the level setpipeline. Moreover, we suggest a new topological EC biomarker which has not been reported in the literature before,moderately correlated with FAZ area with Pearson correlation coefficient of . , which is the end discovery result ofthe radial gradient pipeline. (a) Control: No disease (b) NoDR: Diabetes, but no retinopathy. (c) DR: Diabetes with retinopathy. Figure 13: OCTA scans from Control, DR and NoDR patients. Changes to the microvasculature are apparent withdisease progression. For example, the vessel density reduces and the foveal avascular zone (FAZ), which is the blackregions in the middles of the image, is enlarged and distorted with less circular shape.
We provide the first study of the global topological structure of the OCTA retinal images via the ECC. As mentioenedabove, we take the filtering function to be the pixel intensity value. We calculate the EC for all the possible pixelthreshold values, which in the case of gray-scale images are 256. The ECC is then used as a feature vector forclassification.
W e use for comparison a baseline of biomaker measurements for two knownbiomarkers associated with disease progression, vessel density (VD) and FAZ area, similar to the study in [45], wherethey use 3 biomarkers (VD, FAZ area and vessel calibre). Moreover, we further compare it with the state-of-the-artdeep learning approaches to patient classification. A VGG16 architecture with transfer learning was used as describedin [46] to classify the same OCTA images.
NHS Lothian, Control vs. DRBaseline Our approach VGG16
Overall Acc . ± .
03 0 . ± . ± Sen (Control) ± ± . ± . Spe (Control) . ± .
06 0 . ± . ± AUC . ± . ± ± . ± . ± . ± . . ± . ± ± ± ± . ± . . ± . ± ± Table 2: Table of classification performances in the Control vs. DR studyIn Table 2 we can see the classification performance in the 2-class task (Control vs. DR). Our approach performsbetter than the baseline in all metrics, but the VGG16 method with data augmentation outperforms for accuracy andspecificity. The AUC in both cases is comparable and it has higher variance for VGG16. For OCTAGON, we observethat we achieve a 5% accuracy improvement over the baseline and 3% improvement over the transfer learning approachfor OCTAGON. For the other metrics, the ECC approach is still better than the baseline, and the results are at worstcomparable with VGG16. 17
PREPRINT - F
EBRUARY
17, 2021
We compared the results with the transfer learning approachapplied to the same dataset with data augmentation. The classification statistics with data augmentation for VGG 16 aredisplayed in Table 3.
Controls NoDR DROur approach VGG16 Our approach VGG16 Our approach VGG16ACC ± ± ± . ± . ± ± . ± . ± ± . ± . ± ± ± ± ± ± ± . ± . AUC ± ± ± ± ± . ± . Table 3: Table of classification performances with transfer learningOur results are comparable to the results with data augmentation. Our approach demonstrates good performance for theNoDR class with the highest accuracy of . , followed by DR and then Control. This result signifies its suitability forearly detection tool, as the NoDR images have the slightest of changes compared to DR. The predominant approaches for image analysis in the OCTA literature have centered around a small number ofexplainable candidate biomarkers as suggested by the clinical knowledge acquired [47]. Indeed, quantifiable featurescan be extracted from the OCTA images which are important biomarkers for DR. Statistical studies have identified theusefulness of local and global metrics based on the morphology of the foveal avascular zone (FAZ) and vascular-basedmetrics as biomarkers for distinguishing between healthy and DR eyes. Examples of the former include FAZ area, FAZcontour irregularity [48, 49], while examples of the latter are vessel caliber (VC), fractal dimension (FD), tortuosity,vessel density (VD) and geometric features of the vascular network [41, 50, 51, 52, 53]. For a more detailed review,please refer to [40]. However, there can be up to 25% differences in the measurements of one of the biomarkers knownto be linked to DR (vessel density) and 24% in foveal avascular zone (FAZ) area [54], which is also an early biomarkerfor DR and is enlarged for Diabetic patients.
Vessel Density (VD).
VD is the ratio of the parts of image which are taken by blood vessels to the entire image. VDmeasurements were obtained as described in [55]. For OCTAGON, OOF filter was used, while for NHS Lothian aU-Net approach was adopted due to the availability of manually labelled data.
FAZ area.
The FAZ area is measured by segmenting the FAZ (the black region in the middle as seen in Figure 13a) andcalculating the total area of the segmented region. As there exist different methods for segmenting the FAZ area anddepending on the dataset and the availability of manually segmented data available, the FAZ area is calculated usingtwo different methods. The first one is used for OCTAGON and uses the FAZ segmentation and area calculation asdescribed in [39]. The second one is used for NHS Lothian and follows the methodology as in [55].
In order to identifying a particular area in the image to look at, we would develop a suitable second image for the ECSguided by the accumulated biomedical knowledge of biomarkers VD and FAZ. We would like to focus on either:• a pixel interval, which consists of vessels with pixel values in particular range; or• restrict our attention to the FAZ area in the middle of the image.For these two particular purposes, we will use the ECS with a carefully selected second image M . As we have seenpreviously, this representation is richer. Level set pipeline.
The complement image is the image in which each pixel value is subtracted from the maximumpixel value, 256 in the case of a gray-scale image. In Figure 14c we can see an example of the resulting ECS fromtaking a Control image M in Figure 14a and its complement to be M as shown in Figure 14b. Radial gradient pipeline.
In Figure 15c. we can see an example of the resulting ECS from taking a Control image M in Figure 15a and the second image M in Figure 15b is the radial gradient image . In this case, the choice ofimage M is motivated by the idea to capture the FAZ in the OCTA images, as its enlargement is characteristic for theprogression of DR [54]. By selecting a threshold t for the radial gradient image, we consider the disk with radius t and18 PREPRINT - F
EBRUARY
17, 2021 (a) Control image ( M ) (b) Complement of control image ( M ) (c) ECS( M , M ) Figure 14: Example of the ECS with M being the complement imageits intersection with the OCTA image. Thus, it is in theory possible to detect the FAZ, rendering the radial gradientimage a suitable candidate for a second image M . (a) Control image ( M ) (b) Radial gradient image ( M ) (c) ECS( M , M ) Figure 15: Example of the ECS with M being the radial gradient image Terrain.
The main insight is that there are areas on the ECS which are highly discriminatory between the two groups(Controls and DR). We compute the normalized terrains using all the images. Recall that this is constructed bycomparing point-wise means and standard deviations of the EC. For OCTAGON, we identify 2 regions of interest inthe level set terrain in Figure 16a - regions A and B. For the radial gradient terrain in Figure 16b there is one region Awhich has a high density of red points. The terrains for NHS Lothian are similar, suggesting the existence of underlyingand fundamental topological features, shared between datasets and well captured by the ECS.In Figures 17 we see some examples of the images from the control and DR from Region B, i.e. the the thresholdscorresponding to a point in Region B. The difference between the control and DR is visually clear which shows how theEuler surfaces can yield insights into the differences between classes. In Figure 18, we show images for the sublevel setand radial filtration function. Here the difference is less obvious, however, the difference indicates that the shape of theFAZ (the circular region in the center which is devoid of blood vessels) becomes less circular and enlarged for DR.
Correlate terrain regions with biomarkers.
To establish a link between the topological regions and the knownbiomarkers, we calculate the correlation between the EC of region B in the level set terrain in Figure 16a and vesseldensity (VD) and the correlation between the EC in region A in the radial gradient terrain in Figure 16b and FAZ area.
EC ( p -value), NHS Lothian EC ( p -value), OCTAGON(VD,EC(VD)) . . × − ) 0 . . (FAZ,EC(FAZ)) . . × − ) 0 . . × − ) Table 4: Correlation results for NHS Lothian and OCTAGON
Correlation results.
We can see a summary of the correlation results in Table 4. For FAZ area, the correlation resultsfor NHS Lothian are consistent with the results for OCTAGON. Therefore, there is preliminary evidence that themethod used is robust between datasets from a different device and that correlation between EC(FAZ) and FAZ area19
PREPRINT - F
EBRUARY
17, 2021 (a) Level set terrain with regions A and B (b) Radial gradient terrain with region of interest A
Figure 16: Examples of regions of interest in the terrains of OCTAGONis maintained at the similar levels across devices. However, this is not the case for VD. For OCTAGON there is notsignificant correlation. This could be attributed to the fact that NHS Lothian on average has higher vessel valid visibilityand less motion artifacts [43]. (a) (b) (c) (d)
Figure 17: Examples of level-sets of images of the control (a,b) and DR (c,d) corresponding to Region B. (a) (b) (c) (d)
Figure 18: Examples of sublevel sets with the radial mask shown in yellow with the control shown in (a,b) and DR (c,d)corresponding to Region A. 20
PREPRINT - F
EBRUARY
17, 2021
The Euler characteristic is perhaps the most ubiquitous topological invariant and is certainly one of the easiest tocompute efficiently. Despite its simple nature, it often captures many interesting features of dataset, particularlywhen randomness is involved, where it is one of the few cases where closed-form expressions exist. While the Eulercharacteristic curve has been used extensively in applications for some time now, multi-parameter analogues have notbeen explored. Here we began this exploration, as the general interest in multi-parameter persistence theory shows, theneed to study higher-dimensional parameter spaces is becoming increasingly.Here we showed that considering these Euler surfaces are useful as features for classification tasks, but also for thekind of qualitative analysis that is one of the benefits of TDA. The identification of ranges of parameter values whichdifferentiate models is a very interesting observation for data analysis. These promising results also naturally lead tointeresting mathematical questions which have yet to be explored.1. Can one do interesting analysis on the Euler surfaces. While we have shown the particular results, willa clustering algorithm run on Figure 16, always yield sensible and more importantly meaningful results?Furthermore, are there interesting interpretations of multiple regions of the Euler surface providing gooddifferentiation?2. Our current approach to the Euler surfaces are point-wise, however, can a more functional approach yieldadditional insights. In particular, are there basis functions which accurately capture the “shape" of the curves?In future work, we plan to examine the result of techniques such as functional PCA.3. As mentioned, the Euler surfaces can be thought of as something of a sample or a subset from the Eulercharacteristic transform (ECT). An active area of research is to understand how to evaluate and characterizethe stability of the ECT. While Euler surfaces are not a solution to this – the identification of differentiatingparameter ranges may be one approach to this problem. A completely open question remains how to interpretthis and understand this phenomenon within the broader theory of constructible sheaf transforms.4. We often would like to encode certain invariances into an analysis, e.g. rotational invariance. In the case ofEuler surcase this can be done by the choice of filtering function, as is the case for distance from the center. Anatural extension of this, is can these invariances be extended to the ECT, so that they remain invertible up toan invariance. This seems likely but remains an open question.
References [1] Allen Hatcher.
Algebraic Topology . Cambridge University Press, Cambridge, 2002.[2] Victor Guillemin and Alan Pollack.
Differential topology , volume 370. American Mathematical Soc., 2010.[3] Robert J Adler and Jonathan E Taylor.
Random fields and geometry . Springer Science & Business Media, 2009.[4] Tom Leinster. The euler characteristic of a category. arXiv preprint math/0610260 , 2006.[5] Yuliy Baryshnikov and Robert Ghrist. Target enumeration via euler characteristic integrals.
SIAM Journal onApplied Mathematics , 70(3):825–844, 2009.[6] William D Penny, Karl J Friston, John T Ashburner, Stefan J Kiebel, and Thomas E Nichols.
Statistical parametricmapping: the analysis of functional brain images . Elsevier, 2011.[7] Gunnar Carlsson and Afra Zomorodian. The theory of multidimensional persistence.
Discrete & ComputationalGeometry , 42(1):71–93, 2009.[8] Francesca Cagliari and Claudia Landi. Finiteness of rank invariants of multidimensional persistent homologygroups.
Applied Mathematics Letters , 24(4):516–518, 2011.[9] Martina Scolamiero, Wojciech Chachólski, Anders Lundman, Ryan Ramanujam, and Sebastian Öberg. Multidi-mensional persistence and noise.
Foundations of Computational Mathematics , 17(6):1367–1406, 2017.[10] Masaki Kashiwara and Pierre Schapira. Persistent homology and microlocal sheaf theory.
Journal of Applied andComputational Topology , 2(1):83–113, 2018.[11] Alex McCleary and Amit Patel. Multiparameter persistence diagrams. arXiv preprint arXiv:1905.13220 , 2019.[12] Heather A Harrington, Nina Otter, Hal Schenck, and Ulrike Tillmann. Stratifying multiparameter persistenthomology.
SIAM Journal on Applied Algebra and Geometry , 3(3):439–471, 2019.[13] Matthew Kahle. Topology of random simplicial complexes: a survey.
AMS Contemp. Math , 620:201–222, 2014.21
PREPRINT - F
EBRUARY
17, 2021[14] Omer Bobrowski and Primoz Skraba. Homological percolation and the euler characteristic.
Physical Review E ,101(3):032304, 2020.[15] Skylar Stolte and Ruogu Fang. A survey on medical image analysis in diabetic retinopathy.
Medical imageanalysis , 64:101742, 2020.[16] Keith J Worsley, Jonathan E Taylor, Francesco Tomaiuolo, and Jason Lerch. Unified univariate and multivariaterandom field theory.
Neuroimage , 23:S189–S195, 2004.[17] Rien Van De Weygaert, Gert Vegter, Herbert Edelsbrunner, Bernard JT Jones, Pratyush Pranav, Changbom Park,Wojciech A Hellwing, Bob Eldering, Nico Kruithof, EGP Patrick Bos, et al. Alpha, betti and the megaparsecuniverse: on the topology of the cosmic web. In
Transactions on computational science XIV , pages 60–101.Springer, 2011.[18] Eitan Richardson and Michael Werman. Efficient classification using euler characteristic.
Pattern RecognitionLetters , 49:99–106, 2014.[19] Theresa Heiss and Hubert Wagner. Streaming algorithm for euler characteristic curves of multidimensional images. arXiv: 1705.02045 , 2018.[20] Robert J Adler, Eliran Subag, Jonathan E Taylor, et al. Rotation and scale space random fields and the gaussiankinematic formula.
The Annals of Statistics , 40(6):2910–2942, 2012.[21] Katharine Turner, Sayan Mukherjee, and Doug M. Boyer. Persistent homology transform for modeling shapes andsurfaces.
Information and Inference: A Journal of the IMA , 3(4):310–344, 12 2014.[22] Justin Curry, Sayan Mukherjee, and Katharine Turner. How many directions determine a shape and othersufficiency results for two topological transforms. arXiv preprint arXiv:1805.09782 , 2018.[23] Robert Ghrist, Rachel Levanger, and Huy Mai. Persistent homology and euler integral transforms.
Journal ofApplied and Computational Topology , 2(1-2):55–60, 2018.[24] Yuliy Baryshnikov, Robert Ghrist, and David Lipsky. Inversion of euler integral transforms with applications tosensor data.
Inverse problems , 27(12):124001, 2011.[25] Robert Ghrist and Michael Robinson. Euler–bessel and euler–fourier transforms.
Inverse problems , 27(12):124006,2011.[26] LOU AUTOR VAN DEN DRIES, Lou Van den Dries, et al.
Tame topology and o-minimal structures , volume 248.Cambridge university press, 1998.[27] Masaki Kashiwara and Pierre Schapira. Integral transforms with exponential kernels and laplace transform.
Journal of the American Mathematical Society , 10(4):939–972, 1997.[28] H. Edelsbrunner and J. Harer.
Computational Topology: An Introduction . American Mathematical Society, 2010.[29] T. Ojala, T. Mäenpää, M. Pietikäinen, J. Viertola, J. Kyllönen, and S. Huovinen. Outex-New Framework forEmpirical Evaluation of Texture Analysis Algorithms. In
Proceedings of the 16th International Conference onPattern Recognition , pages 701–706, 2002.[30] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-Based Learning Applied to DocumentRecognition.
Proceedings of the IEEE , 86(11):2278–2324, 1998.[31] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel,Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. the Journal of machine Learning research , 12:2825–2830, 2011.[32] Roger B Nelsen.
An introduction to copulas . Springer Science & Business Media, 2007.[33] Marius Hofert.
Elements of copula modeling with R . Springer, 2018.[34] Dirk P Kroese and Zdravko I Botev. Spatial process generation. arXiv preprint arXiv:1308.0399 , 2013.[35] Herbert Edelsbrunner and Ernst P Mücke. Three-dimensional alpha shapes.
ACM Transactions on Graphics(TOG) , 13(1):43–72, 1994.[36] Kai Yuan Tey, Kelvin Teo, Anna CS Tan, Kavya Devarajan, Bingyao Tan, Jacqueline Tan, Leopold Schmetterer,and Marcus Ang. Optical coherence tomography angiography in diabetic retinopathy: a review of currentapplications.
Eye and Vision , 6(1):1–10, 2019.[37] Jean Marie Ekoé, Marian Rewers, Rhys Williams, and Paul Zimmet.
The epidemiology of diabetes mellitus . JohnWiley & Sons, 2008.[38] National Health Service. Diabetic retinopathy. Available from: , 2020. [Accessed on 1 August 2020].22
PREPRINT - F
EBRUARY
17, 2021[39] Macarena Díaz, Jorge Novo, Paula Cutrín, Francisco Gómez-Ulla, Manuel G Penedo, and Marcos Ortega.Automatic segmentation of the foveal avascular zone in ophthalmological OCT-A images.
PloS one , 14(2), 2019.[40] Xincheng Yao, Minhaj Nur Alam, David Le, and Devrim Toslak. Quantitative optical coherence tomographyangiography: A review.
Experimental Biology and Medicine , 245(4):301–312, 2020. PMID: 31958986.[41] Ylenia Giarratano, Alisa Pavel, Jie Lian, Rayna Andreeva, Alessandro Fontanella, Rik Sarkar, Laura J Reid,Shareen Forbes, Dan Pugh, Tariq E Farrah, et al. A framework for the discovery of retinal biomarkers in opticalcoherence tomography angiography (octa). In
International Workshop on Ophthalmic Medical Image Analysis ,pages 155–164. Springer, 2020.[42] Yali Jia, Ou Tan, Jason Tokayer, Benjamin Potsaid, Yimin Wang, Jonathan J Liu, Martin F Kraus, HrebeshSubhash, James G Fujimoto, Joachim Hornegger, et al. Split-spectrum amplitude-decorrelation angiography withoptical coherence tomography.
Optics express , 20(4):4710–4725, 2012.[43] Xin-Xin Li, Wei Wu, Hao Zhou, Jun-Jie Deng, Meng-Ya Zhao, Tian-Wei Qian, Chen Yan, Xun Xu, and Su-QinYu. A quantitative comparison of five optical coherence tomography angiography systems in clinical performance.
International journal of ophthalmology , 11(11):1784, 2018.[44] Ian A Thompson, Alia K Durrani, and Shriji Patel. Optical coherence tomography angiography characteristics indiabetic patients without clinical diabetic retinopathy.
Eye , 33(4):648–652, 2019.[45] Harpal Singh Sandhu, Nabila Eladawi, Mohammed Elmogy, Robert Keynton, Omar Helmy, Shlomit Schaal, andAyman El-Baz. Automated diabetic retinopathy detection using optical coherence tomography angiography: apilot study.
British Journal of Ophthalmology , 102(11):1564–1569, 2018.[46] Rayna Andreeva, Alessandro Fontanella, Ylenia Giarratano, and Miguel O Bernabeu. Dr detection using opticalcoherence tomography angiography (octa): A transfer learning approach with robustness analysis. In
InternationalWorkshop on Ophthalmic Medical Image Analysis , pages 11–20. Springer, 2020.[47] Ylenia Giarratano, Alisa Pavel, Jie Lian, Rayna Andreeva, Alessandro Fontanella, Rik Sarkar, Laura Reid, ShareenForbes, Dan Pugh, Tariq E. Farrah, Neeraj Dhaun, Baljean Dhillon, Tom MacGillivray, and Miguel O. Bernabeu.A framework for the discovery of retinal biomarkers in Optical Coherence Tomography Angiography (OCTA) .
MICCAI Workshop on Ophthalmic Medical Image Analysis – OMIA 2020 , 2020.[48] Joobin Khadamy, Kaveh Abri Aghdam, and Khalil Ghasemi Falavarjani. An update on optical coherencetomography angiography in diabetic retinopathy.
Journal of ophthalmic & vision research , 13(4):487, 2018.[49] Noriaki Takase, Miho Nozaki, Aki Kato, Hironori Ozeki, Munenori Yoshida, and Yuichiro Ogura. Enlargement offoveal avascular zone in diabetic eyes evaluated by en face optical coherence tomography angiography.
Retina ,35(11):2377–2383, 2015.[50] Minhaj Alam, Yue Zhang, Jennifer I Lim, Robison VP Chan, Min Yang, and Xincheng Yao. Quantitative opticalcoherence tomography angiography features for objective classification and staging of diabetic retinopathy.
Retina ,40(2):322–332, 2020.[51] David Le, Minhaj Alam, Bernadette A Miao, Jennifer I Lim, and Xincheng Yao. Fully automated geometricfeature analysis in optical coherence tomography angiography for objective classification of diabetic retinopathy.
Biomedical optics express , 10(5):2493–2503, 2019.[52] Minhaj Nur Alam, Taeyoon Son, Devrim Toslak, Jennifer I Lim, and Xincheng Yao. Quantitative artery-veinanalysis in optical coherence tomography angiography of diabetic retinopathy. In
Ophthalmic Technologies XXIX ,volume 10858, page 1085802. International Society for Optics and Photonics, 2019.[53] MB Sasongko, TY Wong, TT Nguyen, CY Cheung, JE Shaw, and JJ Wang. Retinal vascular tortuosity in personswith diabetes and diabetic retinopathy.
Diabetologia , 54(9):2409–2416, 2011.[54] Florentina J Freiberg, Maximilian Pfau, Juliana Wons, Magdalena A Wirth, Matthias D Becker, and StephanMichels. Optical coherence tomography angiography of the foveal avascular zone in diabetic retinopathy.
Graefe’sArchive for Clinical and Experimental Ophthalmology , 254(6):1051–1058, 2016.[55] Ylenia Giarratano, Eleonora Bianchi, Calum Gray, Andrew Morris, Tom MacGillivray, Baljean Dhillon, andMiguel O Bernabeu. Automated segmentation of optical coherence tomography angiography images: benchmarkdata and clinically relevant metrics.