## 1 Introduction

Anatomical trees like blood vessels, dendrites or airways, carry information about the organs they are part of, and if we can meaningfully compare anatomical trees and measurements made along them, then we can learn more about aspects of disease related to the anatomical trees [11, 23]. For example, airway wall thickness is known to be a biomarker for Chronic Obstructive Pulmonary Disease (COPD), and in order to compare airway wall thickness measurements in different patients, a typical approach is to compare average airway wall area percentage measurements for given airway tree generations of particular subtrees [11, 10]. These approaches assume that measurements made in different locations of the lung are comparable on a common scale, which is not always the case [11, 10]. If we can compare tree structures attributed with measurements made along them in a way which respects the structure and geometry of the tree, then we can be more robustly compare measurements whose values are location sensitive. In this paper we present a family of kernels for comparing anatomical trees endowed with vector attributes, and use these to get a more detailed understanding of how COPD correlates with airway structure and geometry.

Related work. Several approaches to statistics on attributed (geometric) trees have recently appeared, and some of them were applied to airway trees [7, 8, 17, 16]. These methods only consider branch length or shape and do not allow for using additional measurements along the airway tree, such as branch radius or airway wall area percentage. Moreover, these methods are computationally expensive [6], or need a set of leaf labels [17, 8], making them less applicable for general trees. Sørensen [20] treats the airway tree as a set of attributed branches which are matched and then compared using a dissimilarity embedding combined with a

-NN classifier. The matching introduces an additional computational cost and makes the approach vulnerable to incorrect matches.

*Kernels* are a family of similarity measures equivalent to inner products between data points implicitly embedded in a Hilbert space. Kernels are typically designed to be computationally fast while discriminative for a given problem, and often give nonlinear similarity measures in the original data space. Using the Hilbert space, many Euclidean data analysis methods are extended to kernels, such as classification [4] or hypothesis testing [9]. Kernels are popular because they give computational speed, modeling flexibility and access to linear data analysis tools for data with nonlinear behavior.

There are kernels available for structured data such as strings [5, 13], trees [22], graphs [21, 19, 3] and point clouds [1]. The current state-of-the-art graph kernel in terms of scalability is the Weisfeiler-Lehman (WL) [19] kernel, which compares graphs by counting isomorphic labeled subtrees of a particular type and “radius” . The WL scales linearly in

and the number of edges, but the scalability depends on algorithmic constructions for finite node label sets. Thus, the WL kernel, like most fast kernels developed in natural language processing and bioinformatics

[13, 22, 5], does not generalize to vector-valued branch attributes.Walk- and path based kernels [21, 1, 3], which reduce to comparing sub-walks or -paths of the graphs, are state-of-the-art among kernels which include continuous-valued graph attributes. Random walk-type kernels [21, 1] suffer from several problems including tottering [15] and high computational cost. The shortest path kernel [3] by default only considers path length, and some of the kernels developed in this paper can be viewed as extensions of the shortest path kernel.

Contributions. We develop a family of kernels which are computationally fast enough to run on large datasets, and can incorporate any vectorial attributes on nodes^{1}^{1}1Our formulation allows both node and edge attributes, as edge attributes are equivalent to node attributes on rooted trees: assign each edge attribute to its child node., e.g., shape or airway wall measurements. Using the kernels in classification and hypothesis testing experiments, we show that classification of COPD can be substantially improved by taking geometry into account. This illustrates, in particular, that airway wall area percentage measurements made at different locations in the airway tree are not comparable on a common scale.

We compare the developed kernels to state-of-the-art methods. We see, in particular, that COPD can also be detected from combinatorial airway tree structure using state-of-the-art kernels on tree structure alone, but we show that these contain no more information than a branch count kernel, as opposed to the geometric tree kernels.

## 2 Geometric trees and geometric tree kernels

Anatomical trees like airways are *geometric trees*: they consist of both combinatorial tree structure and branch geometry (e.g., branch length or shape), where continuous changes in the branch geometry can lead to continuous transitions in the combinatorial tree structure. In addition to its geometric embedding, a geometric tree can be adorned with additional features measured along the tree, e.g., airway branch radius, airway wall thickness, airway wall thickness/branch radius, airway wall area percentage in an airway cross section, etc.

###### Definition 1

A *geometric tree* is a pair where is a combinatorial tree with nodes , root and edges , and is an assignment of (geometric) attributes from a vector space to the nodes of , e.g. position or landmark points. An *attributed geometric tree* is a triple where is a geometric tree and is a map assigning a vector valued attribute to each point .

A common strategy for defining kernels on structured data such as trees, graphs or strings is based on combining kernels on sub-structures such as strings, walks, paths, subtrees or subgraphs [5, 13, 22, 21, 19, 3, 1]. These are all instances of the so-called *R-convolution kernels* by Haussler [12]. We shall use *paths* in trees as building blocks for defining kernels on trees.

Let be a geometric tree. Given vertices there is a unique path from to in the tree, defined by the sequence of visited nodes:

where , is the parent node of , more generally , and is the highest level common ancestor of and in . We call the *node-path* from to in and for each let the *node-rootpath* be the node-path from to the root.

If the geometric node attributes denote embeddings of the edge into the ambient space , a continuous path can be defined, connecting the embedded nodes along the embedded tree . We call the *embedded path* from to in .

Throughout the rest of this section, we shall define different kernels for pairs of trees and , where are attributed geometric trees (including non-attributed geometric trees as a special case with ), . All kernels defined in this section are positive semidefinite, as they are sums of linear and Gaussian kernels composed with known feature maps. This is a necessary condition for a kernel to be equivalent to an inner product in a Hilbert space [2], needed for the analysis methods used in Sec. 3.

### 2.1 Path-based tree kernels

All-pairs path kernels. The all-pairs path kernel is a basic path-based tree kernel. Given two geometric trees, it is defined as

(2) |

where is a kernel defined on paths, and , are paths connecting to and to in and , respectively – for instance, and , or and , as defined above. Note that if the path kernel is a path length kernel, then the all-pairs path kernel is a special case of the shortest path kernel on graphs [3].

The kernel should take large values on paths that are geometrically similar, and small values on paths which are not, giving a measure of the alignment of the two tree-paths and , making an overall assessment of the similarity between the two geometric trees and . The all-pairs path kernel is nice in the sense that it takes every possible choice of paths in the trees into account. It is, however, expensive: The computational cost is , where and is the cost of the path kernel .

Rootpath kernels. The computational complexity can be reduced by only considering rootpaths, giving a rootpath kernel defined as:

(3) |

where is a path kernel as before, and is the path from to the root . This reduces the computational complexity to .

### 2.2 Path kernels

The modeling capabilities and computational complexity of the kernels and depend on the choices of path kernel and path representation .

Landmark point representation of embedded paths. From a shape modeling point of view, equidistantly sampled landmark points give a reasonable representation of a path through the tree. Representing paths by equidistantly sampled landmark points , the path kernel is either a linear or Gaussian kernel:

(4) |

for a scaling parameter which regulates the width of the Gaussian.

Node-path kernels. The landmark point kernels are expensive to compute (see Table 1). In particular, two embedded tree-paths may have large overlapping segments without having a single overlapping equidistantly sampled landmark point, as the distance between landmark points depends on the length of the entire path. Thus, most landmark points will only appear in one path, giving little opportunity for recursive algorithms or dynamic programming that take advantage of repetitive structure. To enable such approaches, we use node-paths.

Assume that and are node-paths in , respectively, as defined above, that is, sequences of consecutive nodes in . We define

(5) |

where is a node kernel. In this paper is either a linear kernel without/with additional attributes :

or a Gaussian kernel with/without attributes

where the Gaussian weight parameters are heuristically set to the inverse dimension of the relevant feature space, i.e.,

and .Now the node-rootpath tree kernel can be rewritten as:

(6) |

where . This can be reformulated as a weighted sum of node kernels, giving substantially improved computational complexity:

###### Proposition

For each , let be the set of vertices at level in . Then

(7) |

where is an -dimensional vector whose coefficient counts the number of descendants of at level in , respectively. The complexity of computing is .

When is a linear kernel or , the kernel can be further decomposed as

(8) |

at total complexity / (without/with attributes ). Here, denotes the Kronecker product.

###### Proof

Eq. (7) follows from the fact that the terms in kernel (6) will be counted once for every pair of descendants of and , respectively, which are at the same level. The descendant vectors for all can be precomputed using dynamical programming at computational cost , since , where is defined as left aligned addition of vectors^{2}^{2}2e.g., .. The cost of computing is thus

where is the cost of computing each node kernel .

To prove (7) without attributes , let and be column vectors and use the Kronecker trick ():

The total complexity is thus . Similar analysis proves the attributed case.

### 2.3 Pointcloud kernels

Anatomical measurements can also be weighed by location using position alone in a pointcloud kernel. The pointcloud kernel does not use the tree structure but treats each edges in the tree as a point and compares all points:

(9) |

where is a kernel on attributed edges. We use a Gaussian edge kernel (GPC):

(10) |

The kernel is designed to weight the contribution to the total kernel of the airway wall area percentage kernel value between edges and by the geometric alignment of the same edges, defined by the geometric kernel value .

Embedded paths | Node-path | |

( landmark points) | ||

All-paths | ||

Root-paths | ||

Attributed all-paths | N/A | |

Attributed root-paths | N/A | |

Attributed linear root-paths | N/A | |

Pointcloud kernel | N/A |

### 2.4 Baseline kernels

The kernels presented in this paper are compared to a set of baseline kernels. Standard airway wall area percentage measurements are often compared by using an average measure over parts of the tree or a vector of average measures in chosen generations. We use two baseline airway wall area percentage kernels:

(11) |

(12) |

where is the average airway wall area percentage averaged over all centerline points in the tree, and is a -dimensional vector of average airway wall area percentages averaged over all centerline points in generations in tree . For these kernels (AAW%, AgAW%), linear versions were also computed (i.e. replaced with ), but the corresponding classification results are not reported as they were consistently weaker than the Gaussian kernels.

Airway segmentation is likely more difficult in diseased as opposed to healthy subjects, as also observed by [20]. In order to check whether the number of detected branches may be a bias in the studied kernels, we compare our kernels to a linear and a Gaussian branchcount kernel (LBC/GBC) defined by

(13) |

The linear kernel LBC is the most natural, since the Hilbert space associated to a linear kernel on is just . However, a linear kernel on -dimensional input cannot be normalized, as (14) produces a kernel matrix with entries , and the GBC kernel is used for comparison in Table 4 to show that the geometric tree kernels are, indeed, measuring something other than branch count.

Several state-of-the-art graph kernels were also used. The random walk kernel [21] did not finish computing within reasonable time. The shortest path kernel [3] was computed with edge number as path length, and the Weisfeiler-Lehman kernel [19] was computed with node degree as node label. Results are reported in Tables 2, 3 and 4.

## 3 Experiments

Analysis was performed on airway trees segmented from CT-scans of subjects from a national lung cancer screening trial. Triangulated mesh representations of the interior and exterior wall surface were found using an optimal surface based approach [18], and centerlines were extracted from the interior surface using front propagation [14]

. As the resulting centerlines are disconnected at bifurcation points, the end points were connected using a shortest path search within an inverted distance map of the interior surface. The airway centerline trees were normalized using person height as an isotropic scaling parameter. Airway wall thickness and airway radius were estimated from the shortest distance from each surface mesh vertex to the centerline. The measurements were grouped and averaged along the centerline by each nearest landmark point.

Out of the participants, were diagnosed with COPD level 1-3 based on spirometry, and were symptom free. The minimal/maximal/average number of branches in an airway tree was , respectively.

Kernel | Linear | Gaussian | average | average | Shortest | Weisfeiler |

root-node- | branchcount | AW % | generation | path | Lehman | |

path | AW % | () | ||||

Comp. time | m s | m s | s | s | m s | m s |

### 3.1 Kernel computation and computational time

The kernels listed in table 4 were implemented in Matlab ^{3}^{3}3Software: http://image.diku.dk/aasa/software.php; published software was used for SP, WL [19]. and computed on a 2.40GHz Intel Core i7-2760QM CPU with 32 GB RAM. Each kernel matrix was normalized to account for difference in tree size:

(14) |

An exception was made for linear kernels between scalars (LBC and AAW%), since normalization such kernels results gives matrix coefficients .

Computation times for the different kernels used in the classification experiments in Section 3.3 on airway trees are shown in Table 4. To demonstrate scalability, some of the kernels were ran on

airways from a longitudinal study of the

participants, see Table 2. The slower kernels were not included.For classification and hypothesis testing, a set of airway trees from distinct subjects was used ( diagnosed with COPD at scan time).

### 3.2 Hypothesis testing: Two-sample test for means

Let denote a set of data objects. Given any positive semidefinite kernel there exists an implicitly defined feature map into a reproducing kernel Hilbert space such that for all [2]. Hypothesis tests can be defined in to check whether two samples are implicitly embedded by into distributions on that have, e.g., the same means [9]. Denote by and the sample means of and in

, respectively; we use as a test statistic the distance

between the sample means and check the null hypothesis using a permutation test. Writing

and , we divide into random partitions of size and , , compute the test statistic for each partition, and compare it with the statistic obtained for the original partition . An approximate-value giving the probability of

and coming from distributions with identical means is now given by . The statistic can be computed from a kernel matrix since distances in can be derived directly from the values of using the binomial formula:Using the test with selected kernels we show that healthy airways and COPD airways do not come from the same distributions (Table 3).

Kernel | Gaussian | Gaussian | Average | Generation- |

pointcloud | branchcount | AW-wall % | average AW-wall % | |

-value | ||||

Kernel | Linear | Linear | Shortest | Weisfeiler |

all-node-path | Root-node-path | path | Lehman | |

-value |

### 3.3 COPD classification experiments

Kernel type | Mean class. | Kernel matrix | Mean class. |

accuracy | computation time | accuracy | |

Rootpath, linear (3), (4) | h m s | ||

Rootpath, Gaussian (3), (4) | h m s | ||

All-node-paths, linear (2), (5) | h s | ||

Root-node-path, linear (3), (5) | m s | ||

Root-node-path, Gaussian (3), (5) | h m s | ||

Root-node-path, linear, (3), (5) | m s | ||

airway wall area % attribute | |||

Pointcloud, Gaussian (9) | h m s | ||

Branchcount, linear (13) | s | N/A | |

Branchcount, Gaussian | N/A | ||

Linear kernel on | s | ||

average airway wall area (11) | |||

Gaussian kernel on average | s | ||

airway wall area %, | |||

generations (12) | |||

Shortest path [3] | m s | ||

Weisfeiler Lehman () [19] | m s |

Based on the kernel matrices corresponding to the kernels described in Sec. 3.1 for a set of

airway trees, classification into COPD/healthy was done using a support vector machine (SVM)

[4]. The SVM slack parameter was trained using cross validation on of the entire dataset, and tested on the remaining . This experiment was repeatedtimes and the mean accuracies along with their standard deviations are reported in Table

4. All kernel matrices were combined with the GBC kernel matrix in order to check whether the kernels were, in fact, detecting something other than branch number.## 4 Discussion

We have constructed a family of kernels that operate on geometric trees, and seen that they give a fast way to compare large sets of trees. We have applied the kernels to hypothesis testing and classification of COPD based on airway tree structure and geometry, along with state-of-the-art methods. We show that there is a connection between COPD and airway wall area percentage, and the COPD detected based on our weighted airway wall area percentage kernels is stronger than what can be found using average airway wall area percentage measurements over different airway tree generations, which is commonly done [11, 10].

Efficient kernels for trees with vector-valued node attributes are difficult to design because algorithmically, similarity of vector-valued attributes is more challenging to efficiently quantify than equality of discrete-valued attributes. Nevertheless, some of the defined kernels for vector-attributed trees are fast enough to be applied to large datasets from clinical trials.

Vector-valued attributes are important from a modeling point of view, as they allow inclusion of geometric information such as branch shape or clinical measurements in the trees. However, there is a tradeoff between computational speed and optimal use of the attributes. The efficient node paths are less robust than the embedded paths in airway segmentations with missing or spurious branches, and we observe a small drop in classification performance in Table 4. Rootpath kernels are introduced to improve computational speed. However, they do introduce a bias towards increased weighting of parts of the tree close to the root, which are contained in more root-paths. Gaussian local kernels perform significantly better than linear ones (Table 4), which is particularly pronounced in the pointcloud kernel. In convolution kernels based on quantification of substructure similarity rather than isomorphic substructure, all the dissimilar substructures are still contributing to the total value of the kernel, and the Gaussian local kernel downscales the effect of dissimilar substructures much more efficiently than the linear kernel. This is particularly pronounced in kernels that use geometric weighting of airway wall measurement comparison. Unfortunately, however, algorithmic constructions like the Kronecker trick (Prop. 2.2) do not work for the Gaussian kernels, which do not scale well to larger datasets.

Using hypothesis tests for kernels we show that the healthy and COPD diagnosed airway trees come from different distributions. Using SVM classification we show that COPD can be detected by kernels that depend on tree geometry, tree geometry attributed with airway wall area percentage measurements, or combinatorial airway tree structure. Another efficient detector of COPD is the number of branches detected in the airway segmentation. It is thus important to clarify that our defined kernels are not just sophisticated ways of counting the detected branches. Combining the GBC kernel with the other kernels improves classification performance of the geometrically informed tree and pointcloud kernels, showing that these kernels must necessarily contain independent information, and the connection between COPD and airway shape is more than differences in detected airway branch numbers. In contrast, graph kernels that only use the tree structure are not significantly improved by combination with the branch count kernel. Future work includes efficient ways of computing all-paths kernels with linear node attributes, efficient kernels for trees with errors in them, as well replacing the Gaussian local kernels with more efficient RBF type kernels.

## Acknowledgements

This research was supported by the Danish Council for Independent Research Technology and Production Sciences; the Lundbeck Foundation; AstraZeneca; The Danish Council for Strategic Research; Netherlands Organisation for Scientific Research; and the DFG project ”Kernels for Large, Labeled Graphs (LaLa)”.

## References

- [1] F.R. Bach. Graph kernels between point clouds. In ICML, pages 25–32, 2008.
- [2] C.M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
- [3] K.M. Borgwardt and H.-P. Kriegel. Shortest-path kernels on graphs. ICDM, 2005.
- [4] C.-C. Chang and C.-J. Lin. LIBSVM: A library for support vector machines. ACM Trans. Int. Syst. and Tech., 2:27:1–27:27, 2011. Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
- [5] C. Cortes, P. Haffner, and M. Mohri. Rational kernels: Theory and algorithms. JMLR, 5:1035–1062, 2004.
- [6] A. Feragen. Complexity of computing distances between geometric trees. In SSPR/SPR, pages 89–97, 2012.
- [7] A. Feragen, S. Hauberg, M. Nielsen, and F. Lauze. Means in spaces of tree-like shapes. In ICCV, 2011.
- [8] A. Feragen, J. Petersen, M. Owen, P. Lo, L.H. Thomsen, M.M.W. Wille, A. Dirksen, and M. de Bruijne. A hierarchical scheme for geodesic anatomical labeling of airway trees. In MICCAI (3), pages 147–155, 2012.
- [9] A. Gretton, K.M. Borgwardt, M.J. Rasch, B. Schölkopf, and A.J. Smola. A kernel two-sample test. JMLR, 13:723–773, 2012.
- [10] M. Hackx, A.A. Bankier, and P.A. Genevois. Chronic Obstructive Pulmonary Disease: CT quantification of airways disease. Radiology, 265(1):34–48, 2012.
- [11] M. Hasegawa, Y. Nasuhara, Y. Onodera, H. Makita, K. Nagai, S.i Fuke, Y. Ito, T. Betsuyaku, and M. Nishimura. Airflow Limitation and Airway Dimensions in Chronic Obstructive Pulmonary Disease. Am. J. Respir. Crit. Care Med., 173(12):1309–1315, 2006.
- [12] D. Haussler. Convolution kernels on discrete structures. Technical report, Department of Computer Science, University of California at Santa Cruz, 1999.
- [13] C.S. Leslie and R. Kuang. Fast kernels for inexact string matching. In COLT, pages 114–128, 2003.
- [14] P. Lo, B. van Ginneken, J.M. Reinhardt, and M. de Bruijne. Extraction of Airways from CT (EXACT’09). In 2. Int. WS. Pulm. Im. Anal., pages 175–189, 2009.
- [15] P. Mahé, N. Ueda, T. Akutsu, J-L Perret, and J-P Vert. Extensions of marginalized graph kernels. In ICML, 2004.
- [16] E. Miller, M. Owen, and J.S. Provan. Averaging metric phylogenetic trees. Preprint, http://arxiv.org/abs/1211.7046, 2012.
- [17] T.M.W. Nye. Principal components analysis in the space of phylogenetic trees. Ann. Statist., 39(5):2716–2739, 2011.
- [18] J. Petersen, M. Nielsen, P. Lo, Z. Saghir, A. Dirksen, and M. de Bruijne. Optimal graph based segmentation using flow lines with application to airway wall segmentation. In IPMI, LNCS, pages 49–60, 2011.
- [19] N. Shervashidze, P. Schweitzer, E.J. van Leeuwen, K. Mehlhorn, and K.M. Borgwardt. Weisfeiler-Lehman graph kernels. JMLR, 12:2539–2561, 2011.
- [20] L. Sørensen, P. Lo, A. Dirksen, J. Petersen, and M. de Bruijne. Dissimilarity-based classification of anatomical tree structures. In IPMI, LNCS, pages 475–485, 2011.
- [21] S.V.N. Vishwanathan, N.N. Schraudolph, R.I. Kondor, and K.M. Borgwardt. Graph kernels. JMLR, 11:1201–1242, 2010.
- [22] S.V.N. Vishwanathan and A.J. Smola. Fast kernels for string and tree matching. In NIPS, pages 569–576, 2002.
- [23] G.R. Washko, T. Dransfield, R.S.J. Estepar, A. Diaz, S. Matsuoka, T. Yamashiro, H. Hatabu, E.K. Silverman, W.C. Bailey, and J.J. Reilly. Airway wall attenuation: a biomarker of airway disease in subjects with COPD. J. Appl. Phys., 107(1):185–191, 2009.

Comments

There are no comments yet.