# Comparative Analysis of Merge Trees using Local Tree Edit Distance

Raghavendra Sridharamurthy, , and Vijay Natarajan, R. Sridharamurthy and V.
Natarajan are with the Department of Computer Science and Automation, Indian
Institute of Science, Bangalore, 560012.
E-mail<EMAIL_ADDRESS>
###### Abstract

Comparative analysis of scalar fields is an important problem with various
applications including feature-directed visualization and feature tracking in
time-varying data. Comparing topological structures that are abstract and
succinct representations of the scalar fields lead to faster and meaningful
comparison. While there are many distance or similarity measures to compare
topological structures in a global context, there are no known measures for
comparing topological structures locally. While the global measures have many
applications, they do not directly lend themselves to fine-grained analysis
across multiple scales. We define a local variant of the tree edit distance
and apply it towards local comparative analysis of merge trees with support
for finer analysis. We also present experimental results on time-varying
scalar fields, 3D cryo-electron microscopy data, and other synthetic data sets
to show the utility of this approach in applications like symmetry detection
and feature tracking.

###### Index Terms:

Merge tree, scalar field, local distance measure, persistence, edit distance,
symmetry detection, feature tracking.

##  1 Introduction

Comparative analysis and visualization of scalar fields is an important
problem with applications to feature detection and tracking, symmetry
detection in scalar fields, and in general to the study of time-varying data.
Topological structures like merge trees (Figure 1) provide an abstract and
combinatorial representation of the scalar field. These representations enable
the analysis methods to focus on topological features of interest. A careful
study of similarities and differences between the topology-based
representations can lead to meaningful comparisons of the underlying scalar
fields. Multiple similarity measures (alternatively comparison measures or
distance measures) have been proposed to compare scalar fields and topological
structures [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. However, most of them
describe methods to compare these structures globally. While global
comparisons help us address a variety of interesting problems such as feature
tracking, detection of periodicity in time-varying data, shape comparison, and
temporal summarization, it cannot be used for finer analysis, specifically of
the local structure or substructures of the corresponding topological
features. For example, consider the split trees rooted at $i$ and $j$ in
Figures 2(a) and 2(b). Global comparison measures convey the overall
dissimilarity, but do not capture the similarity between the regions that map
to the pairs of subtrees rooted at $(i_{10},j_{7})$ , $(i_{8},j_{6})$ and
$(i_{7},j_{5})$ respectively. This type of fine grained or multi-scale
analysis leads to interesting applications and hence there is a need for a
measure that detects similarity both locally and across multiple scales.

###  1.1 Contributions

In this paper, we propose a local tree edit distance based method to compare
substructures of scalar fields across multiple scales. The comparison measure
is an adaptation of the global tree edit distance for merge trees (mted)
introduced by Sridharamurthy et al. [13]. However, it is substantially
different in terms of the definition, properties, approach to its computation,
and applications. This paper makes the following key contributions:

  1. 1.

A novel local tree edit distance (lmted) to compare substructures in scalar
fields.

  2. 2.

A proof that it satisfies metric properties.

  3. 3.

A dynamic programming algorithm to compute the lmted efficiently.

  4. 4.

A notion of truncated persistence to compute costs of matching /
correspondences, which brings in the additional benefit of saving computation
time by reducing the number of comparisons.

  5. 5.

Experiments to demonstrate the practical value of the distance towards
symmetry detection at multiple scales, analysis of the effects of smoothing
and subsampling, a fine grained analysis of topological compression, and
applications to feature tracking.

The mted supports only a few of the above-mentioned applications. Even in
these cases, it is restricted to comparisons on a coarser level or requires a
higher level of user intervention. Feature tracking is not possible with mted
without significant modifications.

(a) 2D scalar field

(b) join tree

(c) split tree

Figure 1: Merge trees. (a) A 2D scalar field together with its critical points
and a set of isocontours. (b,c) A merge tree tracks the connectivity of
sublevel sets (preimage $f^{-1}(-\infty,c]$) or the superlevel sets (preimage
$f^{-1}[c,\infty)$). We consistently use the () color map for the scalar field
and the () color map for representing critical points based on their Morse
index (0: minimum, 1: 1-saddle, 2: 2-saddle, 3: maximum).

(a) Split tree 1

(b) Split tree 2

Figure 2: Scalar fields $f_{1}$ and $f_{2}$ that are not globally similar but
contain locally similar regions. Split trees are overlaid on top of the scalar
fields

###  1.2 Related work

Point-to-point comparisons like rms distance or norms such as $L^{p},1\leq
p<\infty$ may be used to compare scalar fields. But, there is no natural
notion of hierarchy that can be harnessed. Further, these measures are
sensitive to small perturbations in the data.

Existing measures that compare topological structures, and by extension the
underlying scalar fields, are global. They sometimes ignore the hierarchy
present in topological structures such as the contour tree and typically
cannot be used to compare local structures in a meaningful way. Sridharamurthy
et al. [13] summarize previously developed comparison measures and describe in
detail several methods including the bottleneck distance between persistence
diagrams [1], Interleaving distance [3], distance based on branch
decomposition of merge trees [5], edit distance between Reeb graphs [10, 2],
functional distortion distance between Reeb graphs [4], and comparing metric
graphs by persistence distortion [8]. Many of these methods are
computationally intractable or they are setup in such a way that the
computation considers the topological structure as a whole without harnessing
the hierarchy present within.

Simple and practical similarity measures that are not metrics have also been
studied. Saikia et al. [7] propose a measure that compares histograms
constructed based on the merge trees. As in the case of bottleneck distance,
this measure ignores the topological structure during comparison. The
applications do involve local structures but they are query based. The measure
can be computed efficiently and is useful in practice. Saikia and Weinkauf
[11] extend this measure and demonstrate applications to feature tracking in
time-varying data. The measure is used to track a single feature (or a
collection of features) if the feature is specified by the user in a query.

Saikia et al. [6] also introduce the extended branch decomposition graph
(eBDG) that describes a hierarchical representation of all subtrees of a merge
tree and designed an efficient algorithm to compare them. They also present
experimental results on time-varying data. While the representation and
comparison is based on the hierarchy, they demonstrate its use in applications
by explicitly choosing or selecting region(s) of interest rather than
considering the collection of all pairs of subtrees.

Scalar fields may also be compared based on their isosurfaces. Since
theoretically, the number of isosurfaces is infinite, the comparison requires
a selection of a finite number of isosurfaces of interest. Tao et al. [14]
extend the notion of isosurface similarity maps, first conceptualized by
Bruckner and Möller [15], to construct matrices of isosurface similarity maps
(MISM), and use it to explore multivariate time-varying data. This involves
construction of self-similarity maps, temporal similarity maps, and variable
similarity maps followed by temporal clustering and variable grouping.
Finally, paths spanning across these maps are used to guide the visual
comparison. The choice of the isovalues used is crucial but the inclusion
hierarchy followed by isosurfaces is not utilised.

Lukasczyk et al. [16, 17] introduce the concept of nested tracking graphs to
track the entire family of isosurfaces over time while preserving their nested
hierarchy of the isosurfaces. Features are often represented by isosurfaces,
so the method applies to feature tracking. While the method does facilitate
tracking features in all scales and across time, it does not support generic
comparison between features.

Symmetry detection in scalar fields is another important application that
involves the comparison of local substructures to decide if the scalar field
contains repeating patterns. The notion of symmetry has been well studied by
Thomas and Natarajan [18, 19, 20]. While the comparison methods that they
describe work well for symmetry detection, the methods have not been applied
to detect local similarities between different scalar fields in general. In
other words, the symmetry detection problem is a special case of the more
general local similarity detection problem.

Sridharamurthy et al. [13] introduce a tree edit distance between merge trees
(mted). The algorithm to compute the distance processes the trees in a bottom-
up fashion. Computing the global distance involves computing distances between
various pairs of subtrees, which are not necessarily merge trees. In this
work, we show how such a global tree edit distance can be extended towards a
fine-grained comparison of scalar fields. Specifically, the global mted is
restricted to cases where only the tree at the top of the hierarchy is
guaranteed to be a merge tree. In contrast, we facilitate comparison between
all pairs of subtrees of merge trees by ensuring that all comparisons are
between trees that are guaranteed to be merge trees.

##  2 Tree Edit Operations

In this section, we introduce necessary background on edit operations between
rooted trees with a focus on merge trees.

###  2.1 Merge tree

Let $f:\mathbb{X}\longrightarrow\mathbb{R}$ denote a scalar function defined
on a manifold domain $\mathbb{X}$. A _level set_ is the preimage $f^{-1}(c)$
of a real value $c$ (Figure 1(a)). The merge tree of $f$ (Figures 1(b) and
1(c)) tracks the connectivity of sublevel sets ($f^{-1}(-\infty,c]$, _join
tree_) or superlevel sets ($f^{-1}[c,\infty)$, _split tree_) [21]. The split
tree (Figure 3) of a generic scalar function is a simple rooted binary tree.
Nodes of the split tree include critical points of $f$: maxima, saddles, and
the global minimum (root). The maxima have 0 children, the saddles have 2
children, and the global minimum has 1 child. All maxima are paired with
saddles based on the notion of topological persistence [22], except for the
global maximum which is paired with the global minimum. A persistence pair
$(m,s)$ represents a topological feature (connected component of superlevel
set) that is created at a maximum $m$ during a downward sweep of the domain
$\mathbb{X}$ and destroyed at a saddle $s$. The _persistence_ of such a
maximum-saddle pair is defined as the difference in function value at the two
critical points, $pers(m)=pers(s)=f(m)-f(s)$. The _persistence diagram_ is a
plot of the persistence pairs $(f(m),f(s))$ on the plane. The join tree is
defined similarly, its nodes consist of minima, saddles, and the global
maximum (root). Several serial and parallel algorithms are available for fast
computation of merge trees [21, 23, 24], [25].

###  2.2 Tree edit distance

The distance between a pair of trees may be defined by introducing edit
operations that transform one tree into another [26]. Each edit operation has
an associated cost. The tree edit distance is defined as the cost associated
with a sequence of edit operations that transforms one tree into another while
minimizing the total cost. We now introduce the edit operations between
labeled trees, following definitions from Zhang [27]. For a node $i$ in a
rooted tree $T$, $deg(i)$ denotes the number of children of $i$ and
$parent(i)$ is its parent in the tree. The maximum degree of a node in $T$ is
denoted as $deg(T)$. The set of labels of $T$ is denoted by $\Sigma$, and
$\lambda\notin\Sigma$ denotes the null or empty character. An empty tree is
denoted by $\theta$. Since we focus on labeled trees, we use the same notation
to refer to the node and its label.

Edit operations. We consider three edit operations that modify a rooted tree,
one node at a time: _relabel, insert,_ and _delete_ , together with associated
costs.

  1. 1.

relabel: A relabel $i\longrightarrow j$ corresponds to an operation where the
label $i\in\Sigma$ of a node is changed to a label $j\in\Sigma$.

  2. 2.

delete: A delete operation $i\longrightarrow\lambda$ removes a node with label
$i\in\Sigma$ and all children of the node $i$ are made children of
$parent(i)$.

  3. 3.

insert: An insert operation $\lambda\longrightarrow j$ inserts a node with
label $j\in\Sigma$ as a child of another node $i$ by moving all children of
$i$ to children of $j$.

We define a cost function $\gamma$ that assigns a non-negative real number to
each edit operation. It is useful if the cost function $\gamma$ satisfies
metric properties [13].

Edit distance. The distance between two trees $T_{1},T_{2}$ is defined as

$\displaystyle D_{e}(T_{1},T_{2})=\min_{S}\\{\gamma(S)\\},$ (1)

where $S$ is a sequence of edit operations that transforms $T_{1}$ to $T_{2}$.
Zhang et al. [28] showed that $D_{e}$ is equal to the cost of an ancestor
preserving mapping between nodes of $T_{1}$ and $T_{2}$ but computing $D_{e}$
is an NP-complete problem.

###  2.3 Constrained edit distance mappings

Incorporating an additional constraint, namely restricting the mapping of
disjoint subtrees to disjoint subtrees, makes the problem computationally
tractable. This constraint is clearly meaningful in the context of merge trees
because disjoint subtrees correspond to spatially disjoint regions in the
domain and hence spatially disjoint features. We now describe a constrained
edit distance $D_{c}$ and the corresponding mappings [27]. The recursive
definition of $D_{c}$ naturally incorporates the constraint. We refer the
reader to the supplementary material for a few illustrative examples.

We use $i$ and $j$ with or without subscripts to denote both the nodes and
their labels in the trees $T_{1}$ and $T_{2}$, respectively. Let
$i_{1},i_{2},\ldots,i_{n_{i}}$ be the children of $i$ and
$j_{1},j_{2},\ldots,j_{n_{j}}$ be the children of $j$. Let $T[i]$ denote the
subtree rooted at $i$ and $F[i]$ denote the unordered forest obtained by
deleting the node $i$ from $T[i]$. Further, let $\theta$ denote the empty
tree. Then,

$\displaystyle D_{c}(\theta,\theta)$ $\displaystyle=0,$ (2) $\displaystyle
D_{c}(F_{1}[i],\theta)$
$\displaystyle=\sum\limits_{k=1}^{n_{i}}D_{c}(T_{1}[i_{k}],\theta),$ (3)
$\displaystyle D_{c}(T_{1}[i],\theta)$
$\displaystyle=D_{c}(F_{1}[i],\theta)+\gamma(i\longrightarrow\lambda),$ (4)
$\displaystyle D_{c}(\theta,F_{2}[j])$
$\displaystyle=\sum\limits_{k=1}^{n_{j}}D_{c}(\theta,T_{2}[j_{k}]),$ (5)
$\displaystyle D_{c}(\theta,T_{2}[j])$
$\displaystyle=D_{c}(\theta,F_{2}[j])+\gamma(\lambda\longrightarrow j),$ (6)
$\displaystyle\scriptsize D_{c}(T_{1}[i],T_{2}[j])$
$\displaystyle=\scriptsize\min\begin{cases}D_{c}(\theta,T_{2}[j])+\min\limits_{1\leq
t\leq n_{j}}\\{D_{c}(T_{1}[i],T_{2}[j_{t}])-D_{c}(\theta,T_{2}[j_{t}])\\},\\\
D_{c}(T_{1}[i],\theta)+\min\limits_{1\leq s\leq
n_{i}}\\{D_{c}(T_{1}[i_{s}],T_{2}[j])-D_{c}(T_{1}[i_{s}],\theta)\\},\\\
D_{c}(F_{1}[i],F_{2}[j])+\gamma(i\longrightarrow j).\end{cases}$ (7)
$\displaystyle\scriptsize D_{c}(F_{1}[i],F_{2}[j])$
$\displaystyle=\scriptsize\min\begin{cases}D_{c}(\theta,F_{2}[j])+\min\limits_{1\leq
t\leq n_{j}}\\{D_{c}(F_{1}[i],F_{2}[j_{t}])-D_{c}(\theta,F_{2}[j_{t}])\\},\\\
D_{c}(F_{1}[i],\theta)+\min\limits_{1\leq s\leq
n_{i}}\\{D_{c}(F_{1}[i_{s}],F_{2}[j])-D_{c}(F_{1}[i_{s}],\theta)\\},\\\
\min\limits_{M_{r}(i,j)}\gamma(M_{r}(i,j)).\end{cases}$ (8)

Here, $M_{r}(i,j)$ is the restricted edit distance mapping between $F_{1}[i]$
and $F_{2}[j]$. Nodes within different trees of $F_{1}$ are mapped to nodes
lying in different trees of $F_{2}$. Essentially, if $M_{r}$ maps node $i_{1}$
to $j_{1}$ and node $i_{2}$ to $j_{2}$, then $i_{1}$ and $i_{2}$ belong to a
common tree in $F_{1}[i]$ if and only if $j_{1}$ and $j_{2}$ belong to a
common tree in $F_{2}[j]$. The minimum cost restricted mapping may be computed
by constructing a weighted bipartite graph (Figure 4) in such a way that the
cost of the minimum weight maximum matching $MM(i,j)$ is exactly the same as
the cost of the minimum restricted mapping $M_{r}(i,j)$,

$\displaystyle\min\limits_{M_{r}(i,j)}\gamma(M_{r}(i,j))=\min\limits_{MM(i,j)}\gamma(MM(i,j))$
(9)

###  2.4 mted between merge trees

The constrained tree edit distance (mted) between two merge trees
$T_{1},T_{2}$ is defined as the constrained edit distance

$\displaystyle\textsc{mted}(T_{1},T_{2})=D_{c}(T_{1},T_{2})$ (10)

where the cost of edit operations is governed by one of the two cost models
(the $L^{\infty}$ cost $C_{W}$ or overhang cost $C_{O}$) introduced by
Sridharamurthy et al. [13] instead of the generic cost model defined by Zhang
[27].

Both the $L^{\infty}$ cost $C_{W}$ and the overhang cost $C_{O}$, are
intuitive and incorporate properties of merge trees. For the sake of
illustration, we describe $C_{W}$ here and refer to [13, Section 4.2.2] for
the definition of $C_{O}$.

Consider nodes $p\in T_{1}$ and $q\in T_{2}$. Both $p$ and $q$ appear as
points $(b_{p},d_{p})$ and $(b_{q},d_{q})$ in the respective persistence
diagrams. Here, $b_{p}$ represents the lower of the function values between
$p$ and its persistence pair, also referred to as the birth time of the
topological feature. Similarly, $d_{p}$ represents the higher value, also
referred as the death time. The $C_{W}$ cost of edit operations are defined as

$\displaystyle\gamma(p\longrightarrow q)$
$\displaystyle=\min\begin{cases}\max(|b_{q}-b_{p}|,|d_{q}-d_{p}|),\\\
\frac{(|d_{p}-b_{p}|+|d_{q}-b_{q}|)}{2}\end{cases}$ (11)
$\displaystyle\gamma(p\longrightarrow\lambda)$
$\displaystyle=\frac{|d_{p}-b_{p}|}{2}$ (12)
$\displaystyle\gamma(\lambda\longrightarrow q)$
$\displaystyle=\frac{|d_{q}-b_{q}|}{2}$ (13)

The constrained edit distance $D_{c}$ with the $C_{W}$ cost model can be used
to compute mted using a dynamic programming algorithm.

##  3 Local Tree Edit Distance

We introduce a local version of $D_{c}$, a local tree edit distance (lmted)
between merge trees. We begin with a few necessary definitions. For
illustration, we consider the split trees in Figure 3. These trees correspond
to the scalar fields in Figure 2.

(a) Split tree 1

(b) Split tree 2

Figure 3: Split trees corresponding to the scalar fields shown in Figure 2.
One persistence pair is shown within each split tree (orange link). When
comparing subtrees $T_{1}[i_{8}]$ or $T_{2}[j_{5}]$, nodes $i_{5}$ and $j_{4}$
are unpaired. So, dummy nodes $i^{\prime}_{10}$ and $j^{\prime}_{7}$ are
inserted with $|f_{1}(i_{10})-f_{1}(i^{\prime}_{10})|<\varepsilon$ and
$|f_{2}(j_{7})-f_{2}(j^{\prime}_{7})|<\varepsilon$ to serve as root and as a
pair of the unpaired node.

###  3.1 Truncated persistence

The cost of edit operations in mted [13] is based on the topological
persistence of the node(s). Specifically, their $L^{\infty}$ cost $C_{W}$ is
the $L^{\infty}$ distance between the point pair in the persistence diagram
corresponding to the two nodes (relabel) or the distance between the point in
the persistence diagram and the diagonal (insert / delete). While comparing
subtrees such as $T_{2}[j_{5}]$ (Figure 3(b)), a node $j_{4}$ whose
persistence pair is the global root $j$ contributes a large value to the edit
distance. Using the same cost for $j_{4}$ while comparing all subtrees
containing the node may not be appropriate. For example, subtrees
$T_{2}[j_{5}]$ and $T_{2}[j]$ map to two regions in the domain that are vastly
different in size. To alleviate this inconsistency, we introduce the notion of
dummy nodes and truncated persistence. While comparing the subtree
$T_{2}[j_{5}]$ with other subtrees, we insert a dummy node $j^{\prime}_{7}$
that serves as a root of the subtree $T_{2}[j_{5}]$ and as the pair of
$j_{4}$. The function value at the dummy node differs from the parent of
$j_{5}$ at most by a small value $\varepsilon$.

An unpaired node in a subtree corresponds to a leaf whose persistence pair is
outside the subtree. We use $i_{u}$ and $j_{u}$ to denote unpaired nodes in
$T_{1}[i]$ and $T_{2}[j]$, respectively. Dummy nodes corresponding to
$T_{1}[i]$ and $T_{2}[j]$ are denoted by $i^{\prime}$ and $j^{\prime}$. For an
unpaired node $i_{u}\in T_{1}[i]$ and a truncated root that is represented by
a dummy node $i^{\prime}$, we define _truncated persistence_ as

$\displaystyle tp_{i^{\prime}}(i_{u})=|f_{1}(i_{u})-f_{1}(i^{\prime})|.$ (14)

Note that a leaf node will be unpaired in some subtree and, in some cases, it
can be unpaired in multiple subtrees (for example, $j_{4}$). The subscript
$i^{\prime}$ in the definition specifies the subtree under consideration.

###  3.2 Truncated cost

The cost of the edit operations need to be updated based on the truncated
persistence for unpaired nodes $i_{u}$ and $j_{u}$. Let $\gamma$ denote the
original cost of an edit operation derived from the $L^{\infty}$ cost $C_{W}$
used in mted. We define a new truncated cost $\gamma^{\prime}$ as

$\displaystyle\gamma^{\prime}(i\longrightarrow
j)=\begin{cases}\gamma(i\longrightarrow\lambda),i\neq
i_{u},j=j_{u}|\lambda,\\\ \gamma(\lambda\longrightarrow
j),i=i_{u}|\lambda,j\neq j_{u},\\\ \gamma(i\longrightarrow j),i\neq
i_{u},j\neq j_{u},\\\ 0,i=i_{u}|\lambda,j=j_{u}|\lambda,\end{cases}$ (15)

###  3.3 lmted between merge trees

We now describe a new local tree edit distance that is appropriate for
localized comparison of merge trees, discuss its properties, and present an
algorithm for computing the distance. The local tree edit distance (lmted) for
a pair of trees rooted at $i$ and $j$ is denoted $\textsc{lmted}(i,j)$ and
defined as follows:

$\displaystyle\textsc{lmted}(i,j)=D^{\prime}(i,j)+\Gamma(i_{u}\longrightarrow
j_{u}).$ (16)

Here, $\Gamma(i_{u}\longrightarrow j_{u})$ denotes the relabel cost computed
using the truncated persistence values of $i_{u}$ and $j_{u}$.
$D^{\prime}(i,j)$ is the modified edit distance between the two trees that
excludes the cost between the unpaired nodes from each tree.

Figure 4: Illustrating lmted. To compute lmted between subtrees rooted at $i$
and $j$, we treat the subtrees containing the unpaired nodes $i_{u},j_{u}$,
labeled as $i_{u_{i}},j_{u_{j}}$ and highlighted in orange, differently. For
other nodes, the formulation is same as mted. For $i_{u},j_{u}$, we use
truncated persistence to determine the costs. In the matching required to
compute $M^{\prime}_{r}(i,j)$, we consider truncated persistence to determine
the weights of all edges incident on $i_{u_{i}},j_{u_{j}}$ (highlighted in
orange). $s$ and $t$ are the source and destination nodes of the flow problem
that is equivalent to the matching problem to determine $M^{\prime}_{r}(i,j)$.
The cardinalities of the two sides are made equal by inserting a set of
$n_{j}$ dummy nodes adjacent to $s$ and $n_{i}$ dummy nodes adjacent to $t$.

We now describe the recursive formulation of $D^{\prime}(i,j)$. The key
difference between lmted and mted is the way in which unpaired nodes are
handled. The original formulation (as for $D_{c}$) applies as is for the
paired nodes, not for the unpaired nodes. The unpaired node changes depending
on the level of the subtree, as we move from the leaves towards the root of
the merge tree. As we traverse a merge tree bottom-up considering all possible
subtrees, a node can be unpaired until we reach the level containing its pair,
following which it remains paired until the end. When the node is unpaired, at
every level its contribution changes. We account for the contribution from the
unpaired node to $D^{\prime}(i,j)$ in a final step and so we are able to
retain a recursive formulation. Once the node is paired, its contribution
(equal to its persistence) does not change further. This demands a new
recursive formulation which can handle both the scenarios. The new formulation
$D^{\prime}$ will thus be a modification of $D_{c}$.

From Figure 4 and using similar terminology as before, let
$i_{1},i_{2},\ldots,i_{n_{i}}$ be the children of $i$ and
$j_{1},j_{2},\ldots,j_{n_{j}}$ be the children of $j$. Let $T_{1}[i]$ denote
the subtree rooted at $i$ and $F_{1}[i]$ denote the unordered forest obtained
by deleting the node $i$ from $T_{1}[i]$. Again, $i_{u}$ and $j_{u}$ are
unpaired nodes in $T_{1}[i]$ and $T_{2}[j]$, respectively. Let $i_{u_{i}}$ be
the child lying on the path between $i$ and $i_{u}$ in $T_{1}[i]$ and
$j_{u_{j}}$ be the child lying on the path between $j$ and $j_{u}$ in
$T_{2}[j]$.

Recall that $\theta$ denotes the empty tree and $D_{c}(,)$ denotes mted. Then
$D^{\prime}$ is recursively defined as follows:

$\displaystyle D^{\prime}(\theta,\theta)$ $\displaystyle=0,$ (17)
$\displaystyle D^{\prime}(F_{1}[i],\theta)$
$\displaystyle=\sum\limits_{k=1,k\neq
u_{i}}^{n_{i}}D_{c}(T_{1}[i_{k}],\theta)+D^{\prime}(T_{1}[i_{u_{i}}],\theta),$
(18) $\displaystyle D^{\prime}(T_{1}[i],\theta)$
$\displaystyle=D^{\prime}(F_{1}[i],\theta)+\gamma^{\prime}(i\longrightarrow\lambda),$
(19) $\displaystyle D^{\prime}(\theta,F_{2}[j])$
$\displaystyle=\sum\limits_{k=1,k\neq
u_{j}}^{n_{j}}D_{c}(\theta,T_{2}[j_{k}])+D^{\prime}(\theta,T_{2}[j_{u_{j}}]),$
(20) $\displaystyle D^{\prime}(\theta,T_{2}[j])$
$\displaystyle=D^{\prime}(\theta,F_{2}[j])+\gamma^{\prime}(\lambda\longrightarrow
j),$ (21)

$\displaystyle min_{T_{2}}=\min\begin{cases}\min\limits_{1\leq t\leq
n_{j},t\neq
u_{j}}\\{D_{c}(T_{1}[i],T_{2}[j_{t}])-D_{c}(\theta,T_{2}[j_{t}])\\},\\\
\\{D^{\prime}(T_{1}[i],T_{2}[j_{u_{j}}])-D^{\prime}(\theta,T_{2}[j_{u_{j}}])\\}\end{cases}$

$\displaystyle min_{T_{1}}=\min\begin{cases}\min\limits_{1\leq s\leq
n_{i},t\neq
u_{i}}\\{D_{c}(T_{1}[i_{s}],T_{2}[j])-D_{c}(T_{1}[i_{s}],\theta)\\},\\\
\\{D^{\prime}(T_{1}[i_{u_{i}}],T_{2}[j])-D^{\prime}(T_{1}[i_{u_{i}}],\theta)\\}\end{cases}$

$\displaystyle min_{F_{j}}=\min\begin{cases}\min\limits_{1\leq t\leq
n_{j},t\neq
u_{j}}\\{D_{c}(F_{1}[i],F_{2}[j_{t}])-D_{c}(\theta,F_{2}[j_{t}])\\},\\\
\\{D^{\prime}(F_{1}[i],F_{2}[j_{u_{j}}])-D^{\prime}(\theta,F_{2}[j_{u_{j}}])\\}\end{cases}$

$\displaystyle min_{F_{i}}=\min\begin{cases}\min\limits_{1\leq s\leq
n_{i},t\neq
u_{i}}\\{D_{c}(F_{1}[i_{s}],F_{2}[j])-D_{c}(F_{1}[i_{s}],\theta)\\},\\\
\\{D^{\prime}(F_{1}[i_{u_{i}}],F_{2}[j])-D^{\prime}(F_{1}[i_{u_{i}}],\theta)\\}\end{cases}$
$\displaystyle D^{\prime}(T_{1}[i],T_{2}[j])$
$\displaystyle=\min\begin{cases}D^{\prime}(\theta,T_{2}[j])+min_{T_{2}},\\\
D^{\prime}(T_{1}[i],\theta)+min_{T_{1}},\\\
D^{\prime}(F_{1}[i],F_{2}[j])+\gamma^{\prime}(i\longrightarrow j).\end{cases}$
(22)

$\displaystyle D^{\prime}(F_{1}[i],F_{2}[j])$
$\displaystyle=\min\begin{cases}D^{\prime}(\theta,F_{2}[j])+min_{F_{j}},\\\
D^{\prime}(F_{1}[i],\theta)+min_{F_{i}},\\\
\min\limits_{M^{\prime}_{r}(i,j)}\gamma^{\prime}(M^{\prime}_{r}(i,j)).\end{cases}$
(23)

All terms that involve subtrees containing the unpaired node (orange) are
updated to incorporate $D^{\prime}$, whereas the constrained edit distance
$D_{c}$ appears elsewhere. The bipartite graph formulation that is used to
compute the minimum cost restricted mapping between forests is also updated to
incorporate $D^{\prime}$ and is now denoted as $M^{\prime}_{r}(i,j)$ (Figure
4, bottom).

###  3.4 Cost model and properties

We can employ either of the costs, the $L_{\infty}$ cost $C_{W}$ or the
overhang cost $C_{O}$ introduced for the mted [13, Section $4.2$]. Both costs
are proven to be metrics. If they remain so even with the newly introduced
cost based on the truncated persistence, then by Zhang [27] the lmted
satisfies metric properties. The proofs of non-negativity and symmetry is
straightforward because $\gamma^{\prime}$ together with $\Gamma$ is defined in
terms of $\gamma$, which in turn satisfies both properties because it is a
combination of sum, max, and min of absolute values.

We now prove the triangle inequality for $\gamma^{\prime}$ together with
$\Gamma$. Let $T_{1}[i]$, $T_{2}[j],$ and $T_{3}[k]$ be three subtrees with
unpaired nodes $i_{u},j_{u},k_{u}$. We insert dummy nodes
$i^{\prime},j^{\prime},k^{\prime}$ and construct trees $T_{1}[i^{\prime}]$,
$T_{2}[j^{\prime}],$ and $T_{3}[k^{\prime}]$. The truncated persistence values
$tp_{i^{\prime}}(i_{u}),tp_{j^{\prime}}(j_{u}),tp_{k^{\prime}}(k_{u})$ in
subtrees $T_{1}[i]$, $T_{2}[j],T_{3}[k]$ are respectively equal to the regular
persistence values within trees $T_{1}[i^{\prime}]$,
$T_{2}[j^{\prime}],T_{3}[k^{\prime}]$. For a given triple $i_{1}\in
T_{1}[i],j_{1}\in T_{2}[j],k_{1}\in T_{3}[k]$, we will show that triangle
inequality holds by considering different cases.
Case 1: Nodes $i_{1},j_{1},k_{1}$ are all unpaired or are all paired. When all
the nodes $i_{1},j_{1},k_{1}$ are paired, $\gamma^{\prime}=\gamma$ and hence
triangle inequality holds. If all are unpaired,
$\gamma^{\prime}(i_{1}\longrightarrow
j_{1})=\gamma^{\prime}(j_{1}\longrightarrow
k_{1})=\gamma^{\prime}(i_{1}\longrightarrow k_{1})=0$. Further,
$\Gamma(i_{1}\longrightarrow j_{1})$, $\Gamma(j_{1}\longrightarrow k_{1})$ and
$\Gamma(i_{1}\longrightarrow k_{1})$ are equal to relabel costs ($\gamma$) for
the trees $T_{1}[i^{\prime}]$, $T_{2}[j^{\prime}]$ and $T_{3}[k^{\prime}]$ and
hence triangle inequality holds [13, Section $4.3$].
Case 2: Two nodes are unpaired. The case where $i_{1}$ and $k_{1}$ are
unpaired while $j_{1}$ is paired is impossible because of the constraint that
unpaired nodes are mapped to unpaired nodes and the operation is forced to be
a relabel.
Case 2.1: $i_{1}$ and $j_{1}$ are unpaired. Then the LHS

$\displaystyle\gamma^{\prime}(i_{1}\longrightarrow
j_{1})+\gamma^{\prime}(j_{1}\longrightarrow k_{1})+\Gamma(i_{1}\longrightarrow
j_{1})$ (24) $\displaystyle=0+\gamma(\lambda\longrightarrow
k_{1})+\Gamma(i_{1}\longrightarrow j_{1})$ (25)
$\displaystyle=\gamma(\lambda\longrightarrow
k_{1})+\Gamma(i_{1}\longrightarrow j_{1}),$ (26)

and the RHS

$\displaystyle\gamma^{\prime}(i_{1}\longrightarrow
k_{1})=\gamma(\lambda\longrightarrow k_{1}).$ (27)

Since $\Gamma(i_{1}\longrightarrow j_{1})\geq 0$ always, we have LHS $\geq$
RHS.
Case 2.2: $j_{1}$ and $k_{1}$ are unpaired, then the LHS

$\displaystyle\gamma^{\prime}(i_{1}\longrightarrow
j_{1})+\gamma^{\prime}(j_{1}\longrightarrow k_{1})+\Gamma(j_{1}\longrightarrow
k_{1})$ (28)
$\displaystyle=\gamma(i_{i}\longrightarrow\lambda)+0+\Gamma(j_{1}\longrightarrow
k_{1})$ (29)
$\displaystyle=\gamma(i_{i}\longrightarrow\lambda)+\Gamma(j_{1}\longrightarrow
k_{1}),$ (30)

and the RHS

$\displaystyle\gamma^{\prime}(i_{1}\longrightarrow
k_{1})=\gamma(i_{1}\longrightarrow\lambda).$ (31)

Since $\Gamma(j_{1}\longrightarrow k_{1})\geq 0$ always, we again have LHS
$\geq$ RHS.
Case 3: A single node is unpaired. This is impossible because of the
constraint that unpaired nodes are mapped only to unpaired nodes via a relabel
edit operation.
In all cases, the truncated cost $\gamma^{\prime}$ together with $\Gamma$
satisfies the triangle inequality. So, it follows from Zhang [27] that lmted
is indeed a metric.

##  4 Computing lmted

We propose a modified dynamic programming based algorithm for computing lmted
between merge trees. The use of truncated persistence as cost of the edit
operations implies that lmted can be computed using solutions to non-
overlapping sub-problems for computing $D_{c}$. However, the dynamic
programming formulation needs to be modified because $D^{\prime}$ recursively
depends both on $D^{\prime}$ and $D_{c}$.

###  4.1 Dynamic Programming tables

Consider the subtrees $T[i_{10}]$ and $T[i_{8}]$ in Figure 3. The node $i_{5}$
is unpaired. Within the subtree $T[i_{8}]$, node $i_{5}$ has truncated
persistence $tp_{i^{\prime}_{8}}(i_{5})$. But, its truncated persistence is
equal $tp_{i^{\prime}_{10}}(i_{5})$ when considering the subtree $T[i_{10}]$.
Dynamic programming works by storing the results of sub-problems within a
table so that it can be reused. Entries corresponding to both (in general,
more than two) values of truncated persistence are required for the
computation.

Consider a subtree $T[i]$ with an unpaired node $i_{u}$ and the path from
$i_{u}$ to the global root $r$ as shown in Figure 5. Let $i_{v}$ be the pair
of $i_{u}$, clearly $i_{v}\in ancestor(i)$. While processing the unpaired node
$i_{u}$, we need to distinguish between two cases:

  1. O:

The contribution of $i_{u}$ is measured by persistence as defined in the usual
sense.

  2. M:

The contribution of $i_{u}$ is measured by an appropriate instance of
truncated persistence.

Case O corresponds to all subtrees rooted at $i_{a}\in ancestor(i)$, where
$i_{v}\leq i_{a}\leq r$ in the directed path highlighted in green in Figure 5.
Case M corresponds to trees rooted at $i_{a}$ such that $i_{u}\leq
i_{a}<i_{v}$, which also includes $i$ as highlighted in red in Figure 5.

Figure 5: Illustration of Cases O and M. The portion of the path colored in
red corresponds to case M while the portion of the path colored in green
corresponds to case O.

We wish to design an algorithm that makes a single pass over the two input
merge trees and computes lmted between all pairs of subtrees. In order to
achieve this objective, we propose a modified dynamic programming method that
uses two tables. One table stores the values of $D_{c}$ for sub-problems, as
defined and proposed for mted [27]. We introduce a second table that stores
$D^{\prime}$, partial solutions to the modified edit distance for sub-
problems. Figure 6 and Figure 7 show the dynamic programming table
corresponding to $D_{c}$ and $D^{\prime}$, respectively. Entries in the second
table $D^{\prime}$ are defined as follows:

  1. 1.

To compute the entry $(i,j)$, we need first populate the entries for the
subtrees of $T_{1}[i]$ and $T_{2}[j]$. Say, we need to compute the entry at
$(i_{k},j_{l})$, where $i_{k}\in T_{1}[i]$, $j_{l}\in T_{2}[j]$. If the
subtrees rooted at $i_{k}$ and $j_{l}$ do not contain the unpaired nodes
$i_{u}$ and $j_{u}$, _i.e.,_ $i_{u}\notin T_{1}[i_{k}]$ and $j_{u}\notin
T_{2}[j_{l}]$, then we pick the corresponding entry from $D_{c}$, namely
$D_{c}(T_{1}[i_{k}],T_{2}[j_{l}])$. Figure 8 and Figure 9 shows the entries
required to compute the entry $(i,j)$ in both tables.

  2. 2.

If the subtrees rooted at $i_{k}$ and $j_{l}$ contain the unpaired nodes, we
refer to the modified table $D^{\prime}$, whose entries are computed using
truncated persistence.

Figure 6: Dynamic programming table for $D_{c}$. The only useful entry in this
table is the one corresponding to $(i,j)$. Other entries do not represent a
meaningful distance because the subtrees compared to compute the entries are
not merge trees. Entries in the row and column labeled $-1$ correspond to
comparison with the empty tree $\theta$. Figure 7: The modified DP table
corresponding to $D^{\prime}$, which is used to compute lmted. This table
contains an additional row and column when compared to the table for $D_{c}$.
This additional row and column stores the truncated persistence values for the
unpaired nodes corresponding to the current level. Entries in this table
represent valid (and useful) distance between merge trees, which are subtrees
of the trees rooted at $i$ and $j$, respectively.

Figure 6 shows the original table used to compute $D_{c}$. It does contain
entries corresponding to different pairs of subtrees, but only the global
entry $(i,j)$ corresponds to a distance between two merge trees. Other local
comparisons involve subtrees that are not merge trees. Figure 7 is the
modified table storing values of $D^{\prime}$. Figures 8 and 9 show the
difference between the original and modified tables, and the entries that are
required to calculate $D_{c}(i_{10},j_{7})$ and $D^{\prime}(i_{10},j_{7})$.
The labels correspond to the merge trees from Figure 3. We denote tree
distance entries by $d_{c}(,)$ and $d^{\prime}(,)$, and forest entries by
$f_{c}(,)$ and $f^{\prime}(,)$. The node $i_{5}$ is unpaired in
$T_{1}[i_{10}]$ and $j_{4}$ is unpaired in $T_{2}[j_{7}]$. In case of $D_{c}$,
there is no distinction between paired and unpaired nodes. So,
$D_{c}(i_{10},j_{7})$ depends only on the following entries of the original
table – $d_{c}(i_{10},-1)$, $d_{c}(-1,j_{7})$, $d_{c}(i_{7},j_{7})$,
$d_{c}(i_{8},j_{7})$, $d_{c}(i_{10},j_{5})$, $d_{c}(i_{10},j_{6})$, and
$f_{c}(i_{10},j_{7})$. In case of $D^{\prime}$, for all subtrees involving
unpaired nodes ($T_{1}[i_{8}],T_{2}[j_{5}]$), the corresponding entries are
read from the modified table, while the remaining entries are read from the
table for $D_{c}$. So, the required entries include $d^{\prime}(i_{10},-1)$,
$d^{\prime}(-1,j_{7})$, $d_{c}(i_{7},j_{7})$, $d^{\prime}(i_{8},j_{7})$,
$d^{\prime}(i_{10},j_{5})$, $d_{c}(i_{10},j_{6})$, and
$f^{\prime}(i_{10},j_{7})$. Further, the entries in the additional row and
column are also required. They contain the truncated persistence values for
the unpaired nodes corresponding to the current level.

Figure 8: Entries required to calculate $D_{c}(i,j)$ () and $D^{\prime}(i,j)$
() in the table for $D_{c}$. The required entry is denoted by a green
square.The lower triangle entries correspond to forests and upper triangle
entries to trees. Figure 9: Entries required to calculate $D^{\prime}(i,j)$ ()
in the modified table. The required entry is denoted by a green square.The
lower triangle entries correspond to forests and upper triangle entries to
trees.

###  4.2 Algorithm and Analysis

Zhang described an algorithm for computing the tree edit distance for labeled
unordered trees [27]. It is a dynamic programming based algorithm that follows
from the properties discussed in Section 2.3. We adapt this dynamic
programming based method for computing lmted but compute and maintain two
tables $D_{c}$ and $D^{\prime}$ simultaneously. The algorithms fills entries
in both tables $D_{c}$ and $D^{\prime}$ iteratively. An entry within a table
is computed and filled if entries corresponding to all sub-problems are
already filled in previously. This implicitly corresponds to traversing the
two input trees in a bottom up fashion in tandem. After both tables are
filled, lmted is computed for all pairs of subtrees following the definition
_i.e._ , as a sum of $D^{\prime}$ and $\Gamma$. In the worst case, the
algorithm computes lmted between all pairs of subtrees in
$O(|T_{1}|\times|T_{2}|\times(deg(T_{1})+deg(T_{2}))\times
log_{2}(deg(T_{1})+deg(T_{2})))$ time, similar to mted. Section 2 of the
supplementary material describes the lmted algorithm in detail with
pseudocode. In practice, the computation is restricted to a much smaller set
of entries that need to be filled within both tables. This restricted set of
entries is identified in a preprocessing step as described in the following
section. Also, note that the time is amortized over all pairs of subtrees.

##  5 Refinement and Optimization

Given two merge trees, we observe that in many applications it is unnecessary
to compute lmted between all pairs of subtrees. This section describes
refinements that reduces the number of pairs of subtrees that are considered
for lmted computation. This optimization leads to faster computation times.
This step may be skipped if it is necessary to compare all pairs of subtrees.
We also describe a refinement that ensures that all subtrees considered for
comparison are merge trees.

###  5.1 Ordering subtrees

The dynamic programming algorithm works for any ordering of nodes. The entries
required to compute the distance between a pair of subtrees are computed if
they are not already available from the table. However, we choose to order the
nodes by assigning a priority based on the size of the subtree rooted at the
node and on the number of grid points in the domain mapped to the subtree
(_i.e._ , volume of the corresponding region in the domain). This ordering
facilitates easy identification of similar regions via visual inspection of
the distance matrix because subtrees of similar size appear in the close
vicinity of one another within the matrix.

###  5.2 Comparison refinement

Since lmted is computed between all pairs of subtrees of the two input trees,
the number of comparisons is determined by the size of the two trees. In
general, the scalar fields to be compared are unrelated, defined on different
domains, and have different ranges. So, comparing all pairs of subtrees is
necessary and unavoidable. However, in several applications, it is not
necessary or meaningful to compare all pairs of subtrees. We describe two such
scenarios to motivate a refinement step that reduces the number of
comparisons.

  1. 1.

Symmetry detection: The scalar field is compared with itself. So, we can
discard comparisons between subtrees with vastly different sizes (for example,
$T_{1}[i_{4}]$ and $T_{1}[i_{10}]$ in Figure 3) or a subtree that is contained
within another subtree (like $T_{1}[i_{11}]$ and $T_{1}[i_{8}]$).

  2. 2.

Time-varying data analysis: While analyzing a time-varying scalar field,
assuming a fine enough temporal resolution, we may discard comparisons between
subtrees with vastly different sizes. We may also discard comparisons between
subtrees that map to regions in the domain with significantly different sizes
(area / volume).

We define a set of criteria used to direct the refinement. Each criterion is a
ratio between measures or a statistic computed for a subtree. Let
$T_{1},T_{2}$ be the two input merge trees with nodes
$\\{i_{1},i_{2},\ldots,i_{n_{i}}\\}$ and $\\{j_{1},j_{2},\ldots,j_{n_{j}}\\}$,
and roots $r_{1}$ and $r_{2}$, respectively. Consider a subtree $T_{1}[i_{k}]$
and the mapping to its associated region in the domain. This region is a
connected component of the preimage of the range of scalar values
corresponding to $T_{1}[i_{k}]$. Assuming that the scalar function is defined
over a 3D domain, the volume of this region may be approximated by counting
the number of sample points (vertices or grid points, depending on whether the
domain is represented using a tetrahedral mesh or a cube grid). The aggregate
persistence $P_{i_{k}}$ of the subtree is computed as the sum of persistence
of all persistence pairs contained within the subtree. Given a pair of
subtrees $T_{1}[i_{k}]$ and $T_{2}[j_{l}]$, we use the ratio between

  1. 1.

number of nodes in the subtrees
$\left|T_{1}[i_{k}]\right|/\left|T_{2}[j_{l}]\right|$,

  2. 2.

volume (or area) of the domain that maps to the two subtrees, and

  3. 3.

aggregate persistence of the subtrees $P_{i_{k}}/P_{j_{l}}$

to determine the refinement. Thresholds for the ratios are determined
empirically. For each criterion, we plot the number of pairs of subtrees
against increasing values of the ratio, identify the value of the ratio
corresponding to a sharp decline or the ’knee’ of the curve, and choose this
value of the ratio as threshold. If the plot does not exhibit a clear knee
then we set the threshold to $0.5$. Further, if $T_{1}=T_{2}$ then we discard
all comparisons between $T_{1}[i_{k}]$ and all subtrees contained within
$T_{1}[i_{k}]$.

###  5.3 Subtree refinement

Next, we preprocess the input to ensure that we compare only those subtrees
that are merge trees.

The subtrees that constitute the sub-problems in the dynamic programming
algorithm for mted [13] are not necessarily merge trees. Consider the scenario
in Figure 3. While trees $T_{1}[i]$ and $T_{2}[j]$ rooted at $i$ and $j$ are
merge trees, $T_{1}[i_{8}]$ and $T_{2}[j_{5}]$ are not merge trees as they
contain unpaired nodes, namely $i_{5}$ and $j_{4}$. The entry
$D_{c}(i_{8},j_{5})$ in the DP table (Figure 6) contains a value that is used
to computed the mted between $T_{1}[i]$ and $T_{2}[j]$. But, it is not a
meaningful distance between subtrees $T_{1}[i_{8}]$ and $T_{2}[j_{5}]$ because
the cost model used for operations related to unpaired nodes $i_{5},j_{4}$
depends on their original persistence. From Figure 2, it is clear that the
regions associated to the two subtrees are similar to each other but the value
of $D_{c}(i_{8},j_{5})$ does not reflect this similarity.

We insert a dummy node $i^{\prime}_{10}$ in $T_{1}[i_{8}]$ with
$|f_{1}(i^{\prime}_{10})-f_{1}(i_{10})|<\varepsilon$ for a small
$\varepsilon>0$. In subsequent computations, we consider
$T_{1}[i^{\prime}_{10}]$ as the merge tree corresponding to the subtree
$T_{1}[i_{8}]$. Similarly, $T_{2}[j_{5}]$ contains the unpaired node $j_{4}$.
We insert a dummy node $j^{\prime}_{7}$ and consider $T_{2}[j^{\prime}_{7}]$
as the merge tree corresponding to $T_{2}[j_{5}]$. This conversion is
consistent with the mapping between subtrees and regions of the domain and
hence results in meaningful distances.

##  6 Applications

In this section, we demonstrate the utility of lmted in applications like
symmetry detection, feature tracking, and spatio-temporal exploration of
scientific data. We also describe results of a comprehensive analysis of the
effects of subsampling, smoothing, and topologically controlled compression.
In all cases where a global comparison is meaningful, results based on mted
[13] are taken as a baseline.

###  6.1 Understanding the local tree edit distance

Figure 10: Understanding lmted. Columns of the distance matrix (DM) represent
subtrees rooted at nodes of merge tree $T_{1}$, rows represent subtrees rooted
at nodes of tree $T_{2}$. Nodes are ordered as per the priority described in
Section 5.1. lmted values are shown using a blue-red colormap ($0$ $0.1$)

.

We begin with a simple study to understand lmted, by comparing two scalar
fields shown in Figures 2, whose split trees are shown in Figure 3. Figure 10
shows the distance matrix (DM), entries corresponding to subtree pairs that
are discarded during the refinement step are blank. We observe two blue blocks
of size $4\times 6$ and $2\times 3$ in the DM, confirming that there are two
sets of similar regions and the pair of similar subtrees
$T_{1}[i_{10}],T_{2}[j_{7}]$. Note that the distances between these similar
regions are very small $\leq 0.000093$ in contrast to the larger value of mted
($=0.33$) between the two trees. Such instances of local similarity without
global similarity is common in scientific data. Further, lmted also captures
similarity at different scales.

###  6.2 Symmetry Detection

(a) Volume rendering of the Rubisco RbcL8-RbcX2-8 complex(EMDB-1654)

(b) Distance Matrix (DM)

Figure 11: lmted values in the DM are shown using a blue-red colormap ($0$
$0.1$)

.

Finding symmetric structures in scalar fields is a challenging problem [18,
19, 20]). The mted driven approach [13] extracts a particular set of high
persistent subtrees that are known to be symmetric and compares them to verify
symmetry. We take a different approach where we detect symmetry directly based
on local similarity by comparing the scalar field with itself. We use CryoEM
data from EMDB [29], which contains 3D electron microscopy density data of
macromolecules, subcellular structures, and viruses. We first compute the
simplified merge tree (using a small persistence threshold $<1\%$) and
consider pairs of subtrees after refinements described in Section 5.

We illustrate and analyse the results using the Rubisco RbcL8-RbcX2-8 complex
(EMDB-1654) shown in Figure 11(a). We compute lmted between all pairs of
subtrees of its merge tree after the refinement step. The resulting DM is
shown in Figure 11(b). Blank regions correspond to subtree pairs that are
discarded during the refinement step. Submatrices highlighted in black
correspond to regions in the data that are symmetric. For clarity, we have
shown these submatrices together with the corresponding regions in Figures 12,
13. We observe that lmted detects symmetric regions at different scales.

Any selection of submatrices from Figure 11(b) with a common color corresponds
to a set of symmetric regions, we highlight some of them. In some cases matrix
entries corresponding to symmetric regions may not appear adjacent to each
other as a submatrix. But, it may be possible to visually identify the entries
as belonging to a single cluster. A row/column reordering helps the
identification of these clusters, see Behrisch et al. [30] for details. The
reordering may be restricted to a chosen submatrix to save computation time.
We describe a few additional experimental results together with row/column
reordering in the supplementary material.

(a) DM $124$

(b) DM $97$

(c) volume $97$

(d) volume $124$

(e) volume $129$

(f) volume $139$

(g) volume $143$

Figure 12: Highlighted submatrices from Figure 11(b) and corresponding
regions. 12(d)-12(g) Regions corresponding to submatrices highlighted in
12(a). 12(c) Region corresponding to submatrix shown in 12(b)

(a) DM $10$

(b) volume $10$

(c) DM $18$

(d) volume $18$

(e) DM $22$

(f) volume $22$

(g) DM $44$

(h) volume $44$

(i) DM $56$

(j) volume $56$

(k) DM $77$

(l) volume $77$

Figure 13: Smaller regions of the Rubisco RbcL8-RbcX2-8 complex (EMDB-1654)
corresponding to the highlighted submatrices.

Comparison with previous methods. Thomas and Natarajan [18] process the branch
decomposition of contour trees by building feature descriptors, and use them
to identify similar subtrees. The main limitation of this approach to symmetry
detection is that it is based exclusively on the structure and may fail when
symmetric regions do not manifest as repeating subtrees. For example, if the
field is noisy, subtrees corresponding to noise have high persistence, or when
the field has large flat regions. Their proposed hierarchy descriptor and
similarity measure is a good estimate but not as accurate as examining the
complete hierarchy. It also ignores the geometry of repeating regions leading
to regions with different geometry grouped together and regions with similar
geometry grouped differently. We use grid points mapped to subtrees, as an
easy-to-compute substitute for geometric information. This also helps us to
find symmetry in multiple scales. We use merge trees instead of contour trees
and avoid computation of extremum graphs, geodesic distances, or contour shape
descriptors in contrast to previous methods [19, 20]. lmted computation is
costly compared to the hierarchy descriptor based comparison [18]. We observe
results similar to previous methods based on explicit geometric shape
descriptors [20], but a theoretical guarantee requires further study.

###  6.3 Analysis of subsampling, smoothing, and topology based compression

We analyse the effects of subsampling, smoothing and topology based
compression [31]. While subsampling and smoothing is applied uniformly across
the domain, the effects of compression vary in different parts of the domain.
We showcase how lmted can be used to analyse these effects meaningfully.

Effects of subsampling and smoothing.  Topology changes due to subsampling and
smoothing are not thoroughly quantified. While previous work does present some
analysis based on the mted, it is global and not capable of providing fine-
grained analysis or explain the non-monotonic variation in many cases. We
present a fine-grained analysis using lmted on a scalar field denoted as
$f_{2}$ [13, Section 5.4] and we use the images of the scalar fields and the
DMs from [13, Figure 14] in Figure 14 to illustrate the benefits.

(a) scalar function $f_{2}$

(b) subsampled $f_{2}$

(c) smoothened $f_{2}$

(d) DM for $f_{2}$, original and subsampled

(e) DM for $f_{2}$, original and smoothened

Figure 14: Measuring the effect of subsampling and smoothing (Images sourced
from Figure 14 from [13, Section 5.4]). 14(a) A synthetic function $f_{2}$
sampled over a $300\times 300$ grid. 14(b) $f_{2}$ subsampled down to a
$30\times 30$ grid over $9$ iterations. 14(c) $f_{2}$ smoothed in $9$
iterations. 14(d) DM showing distance between all pairs of subsampled
datasets. 14(e) DMs showing distances between all pairs of smoothed functions.
Row and column indices correspond to the iteration number, $1$ corresponds to
the lowest resolution/extreme smoothing, $10$ corresponds to the original. We
again use a blue-red colormap (low  high). The scales on colormaps for 14(d)
and 14(e) are different

The non-monotonic variation of the distance along a row / column can be due to
multiple factors. While the subsampling and smoothing affects the number of
critical points, and therefore affects the distance, it is not the only
deciding factor. The distance is also affected by (a) type of critical points
inserted / removed, (b) their function values, and (c) changes in persistence
and pairing. We construct the DMs of the mted and lmted for the subsampled
functions. To highlight the utility of lmted, we pick the non-monotonic
entries indexed $(3,4),(3,5),(3,6)$ from Figure 14(d). The trees are
$|T_{3}|=62,|T_{4}|=66,|T_{5}|=62,|T_{6}|=66$. We observe that
$|T_{3}|=|T_{5}|$ but $D_{c}(T_{3},T_{5})>D_{c}(T_{3},T_{4})$ and
$D_{c}(T_{3},T_{5})>D_{c}(T_{3},T_{6})$ even though
$|T_{3}|\neq|T_{4}|,|T_{3}|\neq|T_{6}|$. Thus, size cannot explain the non-
monotonicity. Also, we notice that $T_{3}$ and $T_{5}$ are structurally
similar, all edits are relabels and there is negligible difference in the
function values of the critical points too. The DMs (Figures 15(a),15(b)) show
small changes in the pattern, but the values are similar. The bottom-right
portions of the DMs along with the values are shown in Figures 15(c), and
15(d). The diagonal entries in left portion of Figure 15(c) related to $(3,4)$
shows a gradual increase, while in case of Figure 15(d) related to $(3,5)$ we
observe an upward spike in the last entry. The corresponding entries for mted
in both cases change gradually, even though for $(3,5)$ the increase is
higher. lmted uses truncated persistence for all subtrees and effect of change
in persistence pairings is seen only in the global comparison, causing a jump.
So, the change in distance means that the subsampling has caused a change in
persistence pairing when we go from resolution $4$ to $5$ and $5$ to $6$ but
no such change when we go from $3$ to $4$. Observation of the pairings
confirms this. We also saw that the pairing changes in $4$ to $5$ was reversed
from $5$ to $6$, thus resulting in a lower value in the entry $(3,6)$.

(a) lmted DM for subtree pair $(T_{3},T_{4})$, size $62\times 66$

(b) lmted DM for subtree pair $(T_{3},T_{5})$, size $62\times 62$

(c) Zoomed DM for $(T_{3},T_{4})$

(d) Zoomed DM for $(T_{3},T_{5})$

Figure 15: Measuring the effect of subsampling using lmted. We use a blue-red
colormap for the distances ($0$ $0.1$). Entries that are discarded due to the
refinement are marked with $1.000$.

Fine-grained analysing using lmted.  lmted can be used in conjunction with
mted to quantify the changes caused by subsampling (or smoothing). This is
achieved by computing mted across all resolutions and checking if the
variation is monotonic. If yes, then the subsampling is likely to have caused
changes only in terms of (a) the number of critical points, (b) the function
values of the critical points, and (c) persistence of critical points. If the
variation is non-monotonic with a jump in the last entries of corresponding
lmted, then irrespective of other factors, there are changes in persistence
pairing resulting in changed matching costs and large changes in distance. Due
to the use of truncated persistence in lmted, we can detect such changes as
jumps in distances. While both mted and lmted may be unstable, we observe in
practice that they are more discriminative than bottleneck and Wasserstein
distances.

Figure 16: Topological effects of compression. Yellow stars in the DM for the
subtree pair $T_{0.5},T_{1}$ correspond to regions that remain unchanged.
Distances in the DM are shown using a blue-red colormap ($0$ $0.4$).

Effect of topologically controlled lossy compression.  Soler et al. [31]
describe a method to compress scalar fields that guarantees topology
preservation. The method ensures that the bottleneck distance between the
persistence diagrams of the compressed and uncompressed field is less than a
user specified threshold. Naturally, the method does not consider spatial or
hierarchical structure since it is restricted to the persistence diagrams. We
present here a fine-grained analysis of the effects of compression using
lmted. Soler et al. employ a topological compression followed by zfp. In our
experiments, we use only the former. We begin by computing merge trees for
both the compressed ($T_{c}$) and uncompressed data ($T_{u}$). Since the two
scalar fields are defined over a common domain, we select pairs of subtrees of
$T_{c}$ and $T_{u}$ that correspond to the same region in the domain and order
them based on region size. We compute lmted between these subtree pairs and
note that as we move up the tree hierarchy, the distance remains $0.0$ for
some pairs. The largest among the pairs represent regions that remain
unchanged post compression. Other lmted values follow a staircase pattern,
staying level for a few pairs followed by a jump in value. The jump indicates
that compression has caused a change in the corresponding subtree. Thus we may
identify and isolate regions where compression has no effect in terms of the
function value followed by regions that are affected, and traversing the
hierarchy of the merge tree lends itself to a multi-scale analysis of the
effects of topological compression.

We show results of our analysis applied on AMP-Activated Protein Kinase
(EMDB-1897). We apply topological compression using compression thresholds
$0.5\%,1\%,2\%$, and compute merge trees $T_{0.5},T_{1},T_{2}$ using TTK [32].
To reduce the tree sizes in the experiment, we consider $T_{0.5}$ as the
baseline uncompressed data. We choose regions with $100\%$ overlap and compute
lmted. In Figure 17, we highlight region(s) that remain unchanged for various
thresholds of compression at multiple scales together with a region that is
affected due to compression. Figure 16 shows the DM for the subtree pair
$T_{0.5},T_{1}$, highlighting unchanged regions by a yellow star. We notice
that 19 regions remain unchanged between $T_{0.5}$ and $T_{1}$, and 3 regions
remain unchanged between $T_{0.5}$ and $T_{2}$. A threshold on lmted may be
used to highlight regions that are either affected or remain unaffected for
various compression thresholds.

Figure 17: CryoEM image of AMP-Activated Protein Kinase (EMDB-1897) using
different compression thresholds - $0.5\%,1\%,$ and $2\%$. The region in red
is the largest region that remains unaffected, the region in light red is the
largest region that remains unaffected for compression threshold of $1\%$, and
orange corresponds to $2\%$. Regions in shades of green are symmetric to the
regions in light red and orange but are affected by the compression. The
regions in shades of blue are also affected by compression. The entire protein
is rendered grey and transparent for context.

###  6.4 Spatio-temporal exploration and feature tracking

Figure 18: Visualizing the top $k$ tracks in the 3D von Kármán vortex street
data. The tracks are generated based on lmted and spatial overlaps, and sorted
based on the weights of the tracks. The top tracks capture the temporal
evolution of a set of primary and secondary vortices.

(a) Query, a primary and a secondary vortex (b) Symmetric vortices

Figure 19: Tracking query regions in the 3D von Kármán vortex street data. A
query containing a primary (green) and secondary (red) vortex in time step 399
is selected 19(a). lmted is used to compute regions (light green and orange)
in the same symmetry class as the query 19(b). All regions in the symmetry
class are tracked backward in time (Top row), and forward in time (Bottom
row).

We demonstrate an application of lmted to time-varying scalar fields, in
particular for identifying and tracking features across time. We consider two
scenarios – identify and track all features to provide an overview and an
interactive query-driven mode for feature tracking.

In order to identify and track all interesting topological features, we begin
by computing a sequence of lmted DMs between consecutive timesteps. We apply
the refinements described in Section 5 and compute spatial overlaps between
regions that correspond to the reduced set of subtree pairs. We construct a
track graph whose nodes represent each region and insert an edge between two
nodes in consecutive timesteps if there is a significant overlap between the
corresponding regions. Long paths in the track graph correspond to long-lived
features. We visualize all long-lived features and their evolution over time
including birth, death, split, and merge events. The individual tracks are
also used as a starting point for further analysis. An alternative approach is
to allow the user to specify one or many features within a particular
timestep. We compute regions that are symmetric to the given feature, compute
tracks for each of these regions, and visualize the tracks.

We demonstrate both scenarios using a 3D Bénard-von Kármán vortex street
dataset. The Okubo-Weiss criterion, indicative of high vorticity regions, is
sampled on a regular grid [11]. The scalar field is available on $192\times
64\times 48$ grid with $508$ timesteps. We compute merge trees for all
timesteps and simplify them using a small persistence threshold of $0.8\%$ to
remove noise. We compute lmted on the simplified merge trees after applying
the appropriate refinement steps mentioned in Section 5. The weight of an edge
in the track graph is set equal to the spatial overlap (volume of overlap
normalized by the volume of union) between the corresponding regions. Overlaps
below a $2\%$ threshold are considered negligible and not included into the
track graph.

We process the track graph to enumerate top tracks ordered by either the
length of the track or the sum of weights (high to low). Further, short tracks
(length $<10$) and tracks whose sum of weights is low ($<3.0$) are removed
from consideration. We observe that the first track is a thin region close to
the cylinder obstruction, which remains almost stationary. Other tracks that
appear at the top of the list include the primary and secondary vortices as
identified by [11], see Figure 18. These vortices are represented as
isosurfaces (isovalue 0.1).

In the second scenario, we use a primary and secondary vortex from timestep
$399$ as a query feature, see Figure 19. First, we compute symmetric regions
within the same timestep in order to highlight other primary and secondary
vortices. Next, we compute tracks that contain the query regions and visualize
them. We observe that in the first step lmted can discriminate between the
primary and secondary vortices and, next, it helps efficiently track the
features (vortices) over time. This demonstrates the utility of lmted in the
exploration of time-varying data.

The accompanying video in supplementary material shows (a) the top tracks
corresponding to the primary and secondary vortices, (b) query based
exploration.

There are a few exceptional situations where lmted is unable to discriminate
between primary and secondary vortices. This happens when, say, the chosen
vortex is a secondary vortex, corresponds to a leaf node in the merge tree,
and matches with a leaf node that corresponds to a primary vortex. Further
spatial overlap tests are necessary to identify that the two regions do not
correspond to each other. To summarize, lmted supports the generation of a
good overview visualization and serves as a starting point for feature
detection and tracking. Subsequent interaction and visualization tasks are
often necessary and these tasks may closely depend on application specific
requirements.

##  7 Conclusions

We described a local comparison measure (lmted) between two scalar fields by
comparing subtrees of their merge trees. The comparison measure supports local
and fine-grained analysis and visualization of similarities and differences
between two scalar fields. The measure satisfies metric properties and can be
efficiently computed. We demonstrate its practical utility via applications to
feature tracking, study of topology controlled compression, and symmetry
identification. In future work, we plan to develop a comparative visualization
framework based on the mted and lmted that may be applied to time-varying and
ensemble data.

## Acknowledgments

We thank Tino Weinkauf for providing us with the 3D vortex street dataset with
Okubo-Weiss criterion in full temporal resolution. This work is supported by a
Swarnajayanti Fellowship from the Department of Science and Technology, India
(DST/SJF/ETA-02/2015-16), a scholarship from MHRD, Govt. of India, and a
Mindtree Chair research grant.

## References

  * [1] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer, “Stability of persistence diagrams,” _Discrete and Computational Geometry_ , vol. 37, no. 1, pp. 103–120, 2007.
  * [2] B. Di Fabio and C. Landi, “Stability of Reeb graphs of closed curves,” _Electronic Notes in Theoretical Computer Science_ , vol. 283, pp. 71–76, 2012.
  * [3] D. Morozov, K. Beketayev, and G. Weber, “Interleaving distance between merge trees,” _Discrete and Computational Geometry_ , vol. 49, no. 52, pp. 22–45, 2013.
  * [4] U. Bauer, X. Ge, and Y. Wang, “Measuring distance between Reeb graphs,” in _Proceedings of the 13th annual Symposium on Computational Geometry_. ACM, 2014, pp. 464–474.
  * [5] K. Beketayev, D. Yeliussizov, D. Morozov, G. H. Weber, and B. Hamann, “Measuring the distance between merge trees,” in _TopoinVis III_. Springer, 2014, pp. 151–165.
  * [6] H. Saikia, H. P. Seidel, and T. Weinkauf, “Extended branch decomposition graphs: Structural comparison of scalar data,” _Computer Graphics Forum_ , vol. 33, no. 3, pp. 41–50, 2014.
  * [7] H. Saikia, H.-P. Seidel, and T. Weinkauf, “Fast similarity search in scalar fields using merging histograms,” in _TopoInVis 2015_ , 2015, pp. 1–14.
  * [8] T. Dey, D. Shi, and Y. Wang, “Comparing Graphs via Persistence Distortion,” _31st International Symposium on Computational Geometry (SoCG 2015)_ , pp. 491–506, 2015.
  * [9] V. Narayanan, D. M. Thomas, and V. Natarajan, “Distance between extremum graphs,” in _PacificVis_ , 2015, pp. 263–270.
  * [10] B. Di Fabio and C. Landi, “The edit distance for Reeb graphs of surfaces,” _Discrete and Computational Geometry_ , vol. 55, no. 2, pp. 423–461, 2016\.
  * [11] H. Saikia and T. Weinkauf, “Global feature tracking and similarity estimation in time-dependent scalar fields,” _Computer Graphics Forum_ , vol. 36, no. 3, pp. 1–11, 2017.
  * [12] R. Sridharamurthy, A. Kamakshidasan, and V. Natarajan, “Edit distances for comparing merge trees,” in _IEEE SciVis Posters_ , 2017.
  * [13] R. Sridharamurthy, T. B. Masood, A. Kamakshidasan, and V. Natarajan, “Edit distance between merge trees,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 26, no. 3, pp. 1518–1531, 2020.
  * [14] J. Tao, M. Imre, C. Wang, N. V. Chawla, H. Guo, G. Sever, and S. H. Kim, “Exploring time-varying multivariate volume data using matrix of isosurface similarity maps,” _IEEE transactions on visualization and computer graphics_ , vol. 25, no. 1, pp. 1236–1245, 2018.
  * [15] S. Bruckner and T. Möller, “Isosurface similarity maps,” _Computer Graphics Forum_ , vol. 29, no. 3, pp. 773–782, 2010.
  * [16] J. Lukasczyk, G. Weber, R. Maciejewski, C. Garth, and H. Leitte, “Nested tracking graphs,” _Computer Graphics Forum_ , vol. 36, no. 3, pp. 12–22, 2017.
  * [17] J. Lukasczyk, C. Garth, G. H. Weber, T. Biedert, R. Maciejewski, and H. Leitte, “Dynamic nested tracking graphs,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 26, no. 1, pp. 249–258, 2019.
  * [18] D. M. Thomas and V. Natarajan, “Symmetry in scalar field topology,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 17, no. 12, pp. 2035–2044, 2011.
  * [19] ——, “Detecting symmetry in scalar fields using augmented extremum graphs,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 19, no. 12, pp. 2663–2672, 2013.
  * [20] ——, “Multiscale symmetry detection in scalar fields by clustering contours,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 20, no. 12, pp. 2427–2436, 2014.
  * [21] H. Carr, J. Snoeyink, and U. Axen, “Computing contour trees in all dimensions,” _Computational Geometry_ , vol. 24, no. 2, pp. 75–94, 2003\.
  * [22] H. Edelsbrunner, D. Letscher, and A. Zomorodian, “Topological persistence and simplification,” in _Foundations of Computer Science_. IEEE, 2000, pp. 454–463.
  * [23] A. Acharya and V. Natarajan, “A parallel and memory efficient algorithm for constructing the contour tree,” in _PacificVis_ , 2015, pp. 271–278.
  * [24] D. Morozov and G. Weber, “Distributed merge trees,” in _Proc. ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , ser. PPoPP ’13, 2013, pp. 93–102.
  * [25] C. Gueunet, P. Fortin, J. Jomier, and J. Tierny, “Task-based augmented merge trees with fibonacci heaps,” in _2017 IEEE 7th Symposium on Large Data Analysis and Visualization (LDAV)_. IEEE, 2017, pp. 6–15.
  * [26] P. Bille, “A survey on tree edit distance and related problems,” _Theoretical Computer Science_ , vol. 337, no. 1-3, pp. 217–239, 2005.
  * [27] K. Zhang, “A Constrained Edit Distance Between Unordered Labeled Trees,” _Algorithmica_ , vol. 15, pp. 205–222, 1996.
  * [28] K. Zhang, R. Statman, and D. Shasha, “On the editing distance between unordered labeled trees,” _Information Processing Letters_ , vol. 42, no. 3, pp. 133–139, 1992.
  * [29] “Protein data bank in Europe,” https://www.ebi.ac.uk/pdbe/emdb/, 2021, accessed: 20-03-2021.
  * [30] M. Behrisch, B. Bach, N. Henry Riche, T. Schreck, and J.-D. Fekete, “Matrix reordering methods for table and network visualization,” _Computer Graphics Forum_ , vol. 35, no. 3, pp. 693–716, 2016.
  * [31] M. Soler, M. Plainchault, B. Conche, and J. Tierny, “Topologically Controlled Lossy Compression,” _IEEE Pacific Visualization Symposium_ , pp. 46–55, 2018.
  * [32] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux, “The Topology ToolKit,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 24, no. 1, pp. 832–842, 2018, https://topology-tool-kit.github.io/.

|  Raghavendra Sridharamurthy is a PhD candidate in computer science at Indian
Institute of Science, Bangalore. He received BE degree in information
technology from National Institute of Technology Karnataka, Surathkal and M.Sc
degree in computer science from Indian Institute of Science. His research
interests include scientific visualization, computational topology and its
applications.
---|---
|  Vijay Natarajan is the Mindtree Chair Professor in the Department of
Computer Science and Automation at Indian Institute of Science, Bangalore. He
received the Ph.D. degree in computer science from Duke University. His
research interests include scientific visualization, computational topology,
and computational geometry. In current work, he is developing topological
methods for time-varying and multi-field data visualization, and studying
applications in biology, material science, and climate science.
---|---