Heuristic Search for Homology Localization Problem and Its Application in Cardiac Trabeculae Reconstruction

1 Introduction

Cardiac trabeculae are fine muscle columns, whose ends attache the ventricular walls (Figure 1(a) [Edwin P. Ewing, 2016]). Those fine structures make up a large portion of left ventricle (23% volume at end-diastolic state), and thus play an important role in diagnosis of cardiac diseases and understanding of cardiac functionality. Modern imaging techniques (e.g., Computed Tomography (CT)) make it possible to capture the trabeculae, within cardiac ventricles (Figure 1(b)). However, it is challenging to accurately extract these complex structures due to the heterogeneity of their intensity and geometry. Existing segmentation models, which typically rely on smoothness and global shape priors [Zhu et al., 1995; Boykov et al., 2001; Cootes et al., 1995], are insufficient for handling this task and tend to remove fine geometric details.

To solve this problem, Gao et al. [2013] proposed a topological method, in which salient topological handles (i.e., thickened loops) are detected as trabeculae. Figure 1(c) shows an example of topological handles in a 2D synthetic image. One major weakness with this method is that the detected handles are not ideally represented. A topological handle can be represented by any loop it contains. For example, the handle in Figure 1(c) can be represented by any of the red, green, and cyan cycles. A sub-optimal representative cycle leads to inaccurate geometric description of the handle, and thus impair downstream analysis. Figure 1(d) shows low-quality reconstruction of trabeculae which are not fully stretched as expected. Wu et al. [2017] showed that by finding the optimal representative cycle of each handle, we can obtain the ideal reconstruction results (Figure 1(e)).

However, computing the optimal representative cycle is indeed a very challenging problem. It has been proven that this problem, called homology localization, is NP-hard, and even NP-hard to approximate within constant ratio [Chen and Freedman, 2010a]. Various approximate algorithms [Chen and Freedman, 2010b] have been proposed, but cannot achieve the optimal result. Known exact algorithms, e.g., [Busaryev et al., 2012], are exponential to the Betti number and thus are impractical. To efficiently solve this problem, we formulate it as a search problem. We assign a binary vector, called the homological annotation, to each edge of the underlying graph. In the graph, a loop represents the handle of interest if and only if its edge annotations sum up to a desired vector. This way, the problem of finding the shortest representative cycle of a handle is reduced into finding the shortest loop in a graph with a given sum of edge annotation.

This problem can be solved using an A* searching framework. Wu et al. [2017] proposed a heuristic function for the A* search. It takes the maximum over a set of heuristic functions, each of which focusing on matching one single bit of the given annotation. However, such scheme ignores the relationship between different bits, and may not be a tight approximation of the real cost. In this paper, we design two novel heuristic functions, which merges multiple bits in edge annotations when calculating the heuristic values, thus leading to a tighter estimation. The experimental results demonstrate that the proposed heuristic functions enable the A* algorithm to reduce the searching time by exploring less nodes, and consume less memory than previous methods.

Our improved hierarchical-clustering heuristic is an example of the disjoint pattern database technique developed in [Korf and Felner, 2002]. The optimal homology cycle has to satisfy the constraint that the sum of all bit-vector representations of the edges in the cycle satisfy a particular homology cycle class. The naive heuristic in [Wu et al., 2017] relaxes the constraint so that each bit sums to its desired value independently. Our improved heuristic clusters the bits into exhaustive and exclusive groups so that each group of bits sums to desired values, thus producing a tighter lower bound than the naive heuristic. Also, the naive heuristic and the optimal homology cycle can be considered two extreme cases of our heuristic, with the group size equal to one or the total number of bits respectively.

The rest of this paper is organized as follows. Section 2 introduces the background theories in algebraic topology. Section 3 formulate the problem and explore its solutions via annotation and covering graph. Section 4 presents the search method and the proposed heuristic function. Section 5 provides the experimental results and discussions.

2 Background

In this section, we introduce the theoretical background of algebraic topology [Munkres, 1984; Edelsbrunner and Harer, 2010]. We are particularly interested in 1-dimensional topological structures, i.e., handles. A d-dimensional simplex or d-simplex σ, is the convex hull of d + 1 affinely independent points. A vertex, edge, triangle and tetrahedron are 0-simplex, 1-simplex, 2-simplex and 3-simplex in topology respectively. A face of a d-simplex is the convex hull of a nonempty subset of its d + 1 vertices. A tetrahedron has four 2-dimensional faces, corresponding to four triangles enclosing the tetrahedron. A triangle has three 1-dimensional faces, i.e., the three edges bounding it. The two end vertices of a edge are its faces. A simplicial complex K = {σ1, σ2, ..., σn} is a collection of simplices satisfying two conditions: 1. any face of a simplex in K is also in K. 2. the intersection of any two simplices in K is either empty or their common face. The dimension of a simplicial complex is the maximum dimension of its elements.

A d-chain c is a formal sum of d-simplices in K, c = ∑σ∈K aσσ, aσ ∈ Z2. Note that Z2 is the binary field with only 0 and 1 as its elements and uses Z2 arithmetic as the addition rule (mod 2). A d-chain c is also a subset of the simplicial complex K, C ⊆ K. A d-chain can be represented by a nd long binary vector, in which n0 is the number of dsimplices in K. The i-th entry is 1 if and only if σi ∈ c. All d-chains form the group of d-chain, Cd(K), which is equivalent to a nd-dimensional binary vector space.

The boundary of a d-simplex is the formal sum of its (d - 1)-faces. The boundary of a d-chain, c, is the sum of boundaries of all d-simplices in c, ∂d(c) = ∑σ∈c ∂d(σ). When the chain is represented by a nd-dimensional vector, the boundary operator is equivalent to a nd-1 × nd binary matrix, whose columns are the boundaries of individual dsimplices. See Figure 2 for an example simplicial complex and its boundary matrices. The boundary of c is equals to multiplying the boundary matrix to c, ∂dc. The group of dboundaries is the image of the (d+1)-dimensional boundary matrix, Bd = im(∂d+1).

A d-cycle is a d-chain with zero boundary. The group of d-cycles is the kernel of boundary operator Zd(K) = ker(∂d). In fact, A d-boundary z is a d-cycle, formally, Bd(K) ⊆ Zd(K). A d-cycle z is a non-boundary cycle if z ∈ Zd(K) - Bd(K).

The group of all d-cycles are partitioned into different homology classes. A homology class is an equivalent set of cycles, whose difference is a boundary. Formally, given a cycle z0 belonging to the class h, the class is h = {z | z = z0 + ∂d+1c, c ∈ Cd+1(K)}. The homology class can also be a coset [z0] = z0 + Bd(K). Any two cycles in the same homology class are homologous. Any cycle in the homology class can be the representative cycle, z0. When z0 is a boundary, [z0] is the boundary group Bd(K). All the ddimension homology classes form the d-dimension homology group Hd(K), which is equal to the quotient of Zd(K) over Bd(K), formally, Hd(K) = Zd(K)/Bd(K). As mentioned before, Cd(K) is isomorphic to a vector space (Z2)nd , where nd is the number of d-simplex in K. In fact, Zd(K), Bd(K), and Hd(K) are all vector space (Z2)nd . The dimension of the homology group Hd(K) is called Betti number, βd. For convenient, we focus on 1-dimensional homology and denote g = β1. A set of g cycles, representing a set of linear independent classes, is called a homology cycle basis.

3 Problem

We formulate the localization problem. Given a simplicial complex and a homology class [z0], the localization problem is to find the shortest cycle within the class [z0], formally,

For convenience, we assume all edges have length 1. length(z) is then the number of edges in a cycle z. We use advanced techniques called the homology annotation to convert the localization problem into a shortest path problem.

3.1 Homology Annotations

Let g be the Betti number of dimension 1. A homology annotation is a mapping from edges to binary vectors with length g, A(e) ∈ {0, 1}g , so that for any two cycles, z1 and z2, their annotations are identical if and only if z1 and z2 belong to the same homology class. The annotation of a path or a cycle is equal to the sum of the edge annotations, A(z) = ∑e∈z A(e). The sum is over Z2 field. The edge annotations can be used to identify whether two cycles belong to the same homology class. For example, Figure 3(a)(right) shows an annotation of the complex in Figure 3(a)(left). An edge with no numbers associated has annotation (0, 0). The cycle (e1, e2, e5, e7) and the cycle (e1, e6, e7) belong to the same homology class, as their annotations are both (0, 1).

The annotation can be computed as follows [Busaryev et al., 2012]: Firstly, a spanning tree was created from the input simplicial complex. Any of the rest of edges which do not belong to that spanning tree creates a cycle by connecting two nodes on the tree. Secondly, a homology cycle basis H was calculated from those created cycles. Finally, the annotation of any edge, e, which does not belong to that spanning tree is calculated as a binary vector x, from a linear function Hx = z, where H is the homology cycle basis, and z is a cycle composed by a path on the spanning tree and the edge e. Any edge belonging to the spanning tree has annotation zero. Figure 3(a) shows the construction of an annotation.

With the annotation, we may reformulate the localization problem as finding the shortest cycle with the same annotation as z0:

Intuitively, we are searching for the shortest cycle with a specific annotation. The annotation can be considered as an advanced constraint one has to satisfy.

3.2 Covering Graph

To compute the shortest cycle within a given homology class, a covering graph with N2g many nodes, where N is the number of vertices of the simplicial complex. The covering graph is built on a computed annotation. See Figure 3(b) for an example. In the covering graph, 2g subgraphs are created by copying the spanning tree for 2g - 1 times. In this covering graph, vertices in each subgraph corresponds to a vertex in the original complex associated with an annotation, e.g., u × (0, 0). Edges with non-zero annotation connects vertices cross different subgraphs. An edge e connects two vertices u ×A1 and v ×A2 if and only if A1 +A(e) = A2.

3.3 Dijkstra's Algorithm

Using the covering graph, we have an algorithm for the localization problem: for a vertex v, find a shortest path from v × (0, . . . , 0) at the subgraph with zero annotation to v ×A(z0) at the subgraph with annotation A(z0). This path is essentially the shortest cycle with the desired annotation, A(z0), that passes v. Repeating this procedure over all possible v's will find the shortest cycle with annotation A(z0).

However, this straightforward algorithm is very expensive. In particular, assume the number of vertices in the original complex is N, the covering graph has size 2gN. Solving the shortest path problem for N times will need to run the Dijkstra's algorithm for N times over the covering graph, with a total complexity of O(2gNg logN).

For the remainder of the paper, we solve this problem using A∗ search with specially designed heuristic functions. With good heuristic functions, we do not need to exhaust the whole covering graph in order to find the shortest path as desired.

4 Heuristic Search

In this section, we present an A∗ algorithm to solve the problem. In particular, we design heuristic functions based on the underlying topology and geometry.

In [Wu et al., 2017], A* algorithm was employed to search on the covering graph to find the optimal cycle. In this work, we propose a novel heuristic function which is much tighter approximation of the real cost compared with the previous one. We first introduce the previous heuristic functions introduced in [Wu et al., 2017]. Let ˜ A = A(z0) be the desired annotation. We focus on the search problem where we have a fixed starting point u = u′ × (0, . . . , 0) and end point ω = u′ × ˜ A in the covering graph, and need to find the shortest path between them. Denote by p(u, ω) the path from node u to node ω, and its cost is c(u, ω). Let A(p) represent the annotation for path p. Denote the heuristic estimation from ω to v as h(ω, v). Then the total estimated cost of the target path is c(u, ω) + h(ω, v), which needs to be minimized.

Here h(ω, v) is the maximum of sub-heuristic functions hi(ω, v), each of which is the estimated shortest distance between ω and v. Such estimation should satisfy the constraint that the i-th bit of the annotation A(p(ω, v)) and ˜ A-A(p(u, ω) are identical. Formally,

where P(ω, v) is the set of paths between ω and v, Ai(α) is the value of the i-th bit of annotation of path α, p(u, ω) is the existing optimal path from u to ω.

4.1 The Proposed Heuristic Function

In order to design a tighter heuristic approximation for A* algorithm, we propose to combine two or three annotation bits in a single heuristic function. Given that each annotation bit corresponds to a certain homology class, intuitively it should be more efficient for A* search employing bit combination when some homology classes/handles are tangled. The overall workflow of our method is as follows. We first analyze the input simplicial complex to obtain an optimal homology cycle basis, i.e., a homology cycle basis whole total length is the shortest [Dey et al., 2010]. Then a hierarchical clustering is applied to cycles in this optimal basis. Note that each cycle in the basis corresponds to one annotation bit. By finding a hierarchical clustering of cycles based on their geometrical proximity, we decide which annotation bits should be grouped together. We use single linkage clustering technique to cluster the bases. We set the maximum number of bit combination to 3. The reason is that more bit combinations would lead to an increase of the expense for computing the heuristic values.

The steps of our method are as follows. First, the optimal cycle bases are calculated from the input complex (see Figure 4). Second, we compute the spatial centers of these basis cycles (see Figure 4(b)). Then a tree is created by iteratively merging two basis cycles whose centers are closest, until only one cycle is left. This process builds clusters and progressively creates a merging tree of cycles. Finally, the constructed tree is iteratively cut from the top level to form several sub-trees, until the number of leaves of all sub-trees is less than or equal to 3. The number of leaves in each remaining subtree decides the number of bits combined in the corresponding heuristic function.

After combining the single-bit heuristics based on the hierarchical clustering strategy. The set of g bits of the annotation vector is partitioned into m ≤ g disjoint sets, s1, · · · , sm, so that {1, . . . , g} = s1 ∪ s2 ∪ · · · ∪ sm. The heuristic function can be written as a maximum of m sub-heuristic functions, each of which is constrained by a set of bits, si. Formally, we

When si only contains a single bit, hsi is a single-bit function as in Equation (4.2). Others are combinations of k bits. In experiments, we set k ≤ 3.

Those proposed heuristic functions are able to reduce the searching time, due to theirs intrinsic tighter estimation. In particular, for our problem, the time consumption of A* search mainly consists of three parts, including the time for creating the covering graph, calculating the heuristic value and searching. While the new heuristic function improves search efficiency (thus decreases search time), it is more timeconsuming to compute. In particular, to compute a single-bit heuristic function, hi, as in Equation (4.2) is equivalent to compute the shortest path in a graph of size 2N. The graph is built similar to the covering graph. We use the original graph and a copy of it, corresponding to A = (· · · , 0, · · · ) and A = (· · · , 1, · · · ), respectively. The i-th bit with 0 or 1 is the bit of interest. To run Dijkstra's algorithm on this 2N size graph takes O(N logN) time. In the new proposed method, a particular heuristic function hsi may involve k different bits of interest. To compute it requires building a size 2kN graph, consisting of 2k different copies of the original graph, with different combinations of annotations on the k bits of interest. The computation is then O(2kN(k log 2 + logN)). In other words, our method is able to achieve a better trade-off among different factors and therefore obtain better practical performance. Using smaller k will save heuristic function computation time. Using larger k will lead to tighter bound and thus save search time. We will evaluate the cases when k is upperbounded by 2 and 3, respectively.

The proof of admissibility is relatively straightforward; the problem is essentially a shortest path problem over the covering graph with constraints on the allowable paths (specified in terms of homology annotation). Each heuristic function is the shortest path distance to the destination with only partial constraints (a few bits of the annotation) satisfied. Each heuristic function hsi is naturally a lower-bound of the true shortest path distance. The admissibility will then follow. Also, because the shortest paths over the covering graph are essentially distances in an euclidean space, which automatically satisfy the triangular inequality, the proof of consistency immediately follows. Since the new heuristic combines several bits for heuristic function, it is guaranteed to provide a tighter bound than the previous one (which takes the maximum of all single-bit heuristic functions). We will demonstrate the improvement in experiments.

5 Experimental Results and Discussion

We evaluate the performance of the A* search with proposed heuristic functions. After calculating the persistence diagram, we compare different searching strategies to locate the optimal homology cycles on synthetic and real images, including (1) the NaiveTopo [Gao et al., 2013], which produces non-optimal cycles; (2) OptTopoDij [Wu et al., 2017], which employs Dijkstra's algorithm on searching; (3) OptTopoA* [Wu et al., 2017], which uses A* algorithm with less tighter heuristic function; (4) the proposed methods, which use tighter heuristic functions.

We measure the performance of the following metrics: number of nodes expended (Nodes), and the total time and memory consumption (Mem). In particular, for the total (Total) time cost we further define the following fine-grained metrics: the time on constructing the covering graph (Graph), calculating the heuristic value (H Val) and searching (Search) with A* algorithm, where those terms in the parenthesis, are used from Table 1. Dij, A*, method (1), and method (2), represent the OptTopoDij, OptTopoA*, proposed 2-bit annotation combination in heuristic function, and proposed 3-bit annotation combination in heuristic function, respectively. We conduct the experiments on Intel i5-3427U CPU with 8GB memory.

5.1 Synthetic Experiments

We use randomly generated 2D synthetic images of size 300 × 300 (see Figure 5). The synthetic images have different numbers of cycles, from 3 (easy) to 7 (hard task). Experimental results are shown in Table 1. It can be observed that, the A* based searching methods, including OptTopoA*, and the proposed method (1) and (2), expand less number of nodes than OptTopoDij, and consume less memory for all the synthetic images.

For creating the covering graph, the time consumption of OptTopoDij is more than that of other methods, and is larger when the cycle number is larger (e.g., on harder task). The reason is that the size of covering graph created by OptTopoDij is exponential to the Betti number, g. However, A* based searching methods only create the covering graph with vertices belonging to the edge of specific homology class, and thus consume much less time.

OptTopoDij does not need to calculate heuristic values. The proposed method (2) shows better performance on calculating the heuristic. The reason is method (2) uses the heuristic function with 3-bit combination on annotation. The proposed method (1) and method (2) are only effective when the image has sufficiently complex topology. As a result, the advantage of time consumption of method (2) on calculating the heuristic values, varies depending on the images. The searching time consumption mainly depends on how many nodes explored. A* based methods always have better performance, especially method (2).

We can see the proposed heuristic function drastically improved the performance. For the total time consumption and node expended, the speedup of the proposed methods are 2.8, 1.8, 1.3, 1.5 and 1.6 for the 3-, 4-, 5-, 6- and 7-cycle images, respectively. The number of nodes are reduced by 5.6, 1.2, 2.1, 1.2 and 1.7 times, respectively. Except for the 3-cycle case, the improvement is consistent even for larger problems. Those results are the average performance. If we only focus on the most difficult handle (i.e., the one requiring the most time to optimize), the improvement is even more significant. For the 3-, 4-, 5-, 6- and 7-cycle images, the speedup on the most difficulty handle are 3.4, 2.6, 1.6, 4.7 and 3.8, respectively. The number of expanded nodes are reduced by 5.8, 1.3, 16.1, 7.6 and 2.4 times, respectively.

This selective way of combining bits is critical. We can show that it is much more effective than randomly combining bits. For synthetic input images, we evaluate the baseline of randomly combining bits. For the 3-cycle synthetic image, using the baseline of combining two random bits, the total run time and number of expanded nodes are 73.51 seconds and 47769 respectively. This baseline is much worse than the proposed methods, and are even worse than the previous A* heuristic method. For the 7-cycle synthetic image, the baseline of combining two random bits has 280.52 seconds total time and 141309 expanded notes. It is even worse than the previous A* heuristic method. Randomly combining 3 bits only performs worse.

5.2 Real World Data

We also validate the proposed method on trabeculae segmentation on cardiac CT images (30×17×22), as shown in Figure 6. A standard segmentation method (Figure 6(a)) fails to segment the trabeculae. Although the NaiveTopo (Figure 6(b)) is able to identify the trabeculae, it produces results with wiggling morphology due to the non-optimal homology cycles obtained. In contrast, heuristic search algorithms generate high quality segmentation (Figure 6(c)). Also, the proposed methods achieves better performance in total time consumption and node expanded among different A* based searching methods, as shown in Table 1. Note that for the real world image, Dijkstra's algorithm cannot finish. Only heuristic search methods can find optimal cycles. The experimental results indicate that the A* search with proposed heuristic functions, consume less time and explore less nodes than the existing methods, while identify the same topological handles.

6 Conclusion

We propose a novel heuristic function for the computation of optimal representative cycles of topological handles. The method is validated on synthetic data and applied to real world cardiac image data. Topology-inspired methods can be useful in other applications where we need to reconstruct fine-scale structures with large variation in geometry and appearance. Examples include but are not limited to neurons [Uzunbas et al., 2016; Li et al., 2017] and vessels [Rudyanto et al., 2014].

More advanced topological information such as persistent homology can be used in training classifiers [Chen et al., 2019] and improving clustering [Ni et al., 2017]. In these contexts, optimal persistent homology cycles [Wu et al., 2017] add geometric information to topological descriptors and can improve analytical power.

Acknowledgments

We thank reviewers for their insightful comments and suggestions. We thank Dr. David Mount for suggesting us to look into this problem deeper. This work was partially supported by NSF IIS-1855759, CCF-1855760, CCF-1733843, IIS-1829560, CCF-1733798 and CCF-1740761.