---
title: Distributed Hierarchical Contour Trees
url: https://www.emergentmind.com/topics/distributed-hierarchical-contour-trees
type: topic
---

# Distributed Hierarchical Contour Trees

Searching arXiv for the specified paper and closely related contour-tree work to ground the article in published sources.
Distributed hierarchical contour trees (DHCTs) are distributed representations of contour trees for large-scale scientific data, intended for distributed computation and storage in high-performance computing settings. In the formulation summarized in "Distributed Augmentation, Hypersweeps, and Branch Decomposition of Contour Trees for Scientific Exploration" [2408.04836], a DHCT is not only a storage scheme for the global contour tree but also a substrate for augmentation, hypersweeps, branch decomposition, and contour extraction. The central problem addressed is that, although recent work introduced distributed hierarchical contour trees for distributed computation and storage of contour trees, effective use of these distributed structures in analysis and visualization requires subsequent computation of geometric properties and branch decomposition to support contour extraction and exploration [2408.04836].

## 1. Distributed representation and hierarchy

Let $M$ be a $d$-dimensional mesh of $N$ cells and $N$ vertices. The global contour tree is written as $T=(V,A)$, where $V$ is the vertex set and $A$ is the arc set. The representation used for DHCTs distributes $T$ across $P$ MPI ranks by assigning each rank $r$ a contiguous block of the mesh $M_r \subset M$, with blocks forming a decomposition $M=\bigcup_{r=0}^{P-1} M_r$ and overlap on block boundaries [2408.04836].

On each rank, the stored object is a local sub-tree $T_r=(V_r,A_r)$ such that
$$
V_r=\{v\in V \mid \text{the topological zone (preimage) of } v \text{ intersects } M_r\},
$$
and
$$
A_r=\{a\in A \mid \text{the zone of } a \text{ intersects } M_r\}.
$$
Across all ranks, $\bigcup_r V_r = V$ and $\bigcup_r A_r = A$, so the global contour tree is recovered as the union of the local sub-trees, with consistent global identifiers. A process-assignment function
$$
p: V \cup A \to \{0,\ldots,P-1\}
$$
specifies ownership, with $p(x)=r$ when $x$ is stored on rank $r$ as part of $M_r$'s local sub-tree. On shared arcs or shared nodes, the descriptor is replicated on each rank in the union.

The hierarchy arises from a recursive fan-in of merge operations over $\log P$ rounds, followed by a fan-out reinsertion of suppressed interior branches. Rounds are labeled $\ell=0,1,\ldots,R=\lceil \log_2 P \rceil$. At level $\ell$, ranks are paired or grouped to exchange and merge contour-tree fragments. In the final DHCT, each vertex $v$ stores an integer $\mathrm{level}_\ell(v)$ indicating the level at which it was introduced or reinserted during the fan-out. The vertex set is further decomposed as $V = V_{\mathrm{crit}} \cup V_{\mathrm{reg}}$, where $V_{\mathrm{crit}}$ is the set of critical points and $V_{\mathrm{reg}}$ the set of inserted regular nodes used to subdivide arcs for correct prefix-sum ordering; each arc $a \in A$ is annotated by the level $\ell(a)$ at which it was last modified.

This representation makes the contour tree simultaneously global and rank-local. A plausible implication is that the DHCT hierarchy is designed to preserve enough structural redundancy for parallel topological queries, while still aligning the topology with a mesh decomposition natural to MPI-style execution.

## 2. Augmentation of lazy attachment points

After the DHCT has been constructed through the fan-in/fan-out procedure of Carr et al., some regular attachment points remain lazy: they are recorded only on the rank owning the child subtree and are not inserted on the parent arc across all ranks. The augmentation algorithm inserts each such point as an explicit node along its parent superarc on every rank that holds that superarc [2408.04836].

The algorithm has two phases. In the fan-in exchanges, for levels $\ell=1,\ldots,R$, all attachment points on arcs introduced at level $\ge \ell$ are identified and exchanged between the two partner ranks of that round. After $R$ rounds, every rank holding a given arc knows all of its to-be-inserted attachment points. In the fan-out insertion, for levels $\ell=R,\ldots,1$, the local DHCT is reconstructed by splitting each arc $a$ at level $\ell$ into a sequence of arcs through the sorted list of new attachment points.

The exchanged data are stored as attachment-point lists $AP_r[\ell]$, where $AP_r[\ell]$ is a list of $(\mathrm{arc\_id}, \mathrm{vertex\_id}, f\text{-value})$ for each lazy attachment on arcs with level $\ge \ell$ that originate from rank $r$'s local suppressed branches. At each level, a partner is defined by $r \;\mathrm{XOR}\; 2^{\ell-1}$; lists are exchanged with `MPI_Sendrecv`, merged, and deduplicated. In fan-out insertion, each surviving arc $a$ is replaced by a chain of arcs through the attachment points sorted by $f$-value, after which vertices in the new local DHCT are renumbered consecutively.

Communication cost per level is $O(|AP_r[\ell]|)$ words. Because each attachment originates from a boundary between blocks, the typical three-dimensional case gives $|AP_r[\ell]| = O(N^{2/3}/P^{2/3})$ per round. Over $R=O(\log P)$ rounds, the communication volume is
$$
\sum_{\ell=1}^R O(N^{2/3}/P^{2/3}) = O((N^{2/3}/P^{2/3}) \log P).
$$
Local work per level is dominated by merging and sorting $AP_r[\ell]$, at cost $O(|AP_r[\ell]| \log |AP_r[\ell]|)$. The report notes that, in the worst case, one could have $O(t)$ attachments, giving $O(t \log t \cdot R)$ total work, but that in practice $t=O(N^{2/3})$ and attachment-point lists remain small.

A common misconception is that distributed storage of contour trees is by itself sufficient for downstream scientific exploration. The augmentation procedure directly contradicts that view: distributed prefix operations, called hypersweeps, require those lazy attachment points to be inserted explicitly on every relevant rank.

## 3. Hypersweeps and distributed subtree measures

Once the DHCT is fully augmented, every rank has, for each hyperarc $h$, the exact ordered list of its supernodes. A hyperarc is defined here as a maximal monotone chain of superarcs directed toward the global root. This makes it possible to compute any associative subtree measure by a distributed prefix-sum, termed a hypersweep, plus a fan-in reduction [2408.04836].

The procedure has three steps. First, each rank computes a per-arc measure $m(a)$ on each arc $a \in A_r$. The explicit example given is
$$
m(a)=|Z(a)\cap M_r|,
$$
the number of regular mesh vertices in the topological zone $Z(a)$. Second, per-arc measures are converted to per-hyperarc prefix sums
$$
s(h,i)=\sum_{j \le i} m(a_j)
$$
by an $O(\log |h|)$ shared-memory parallel scan across the $k$ arcs in $h$. Third, a distributed fan-in reduction accumulates each hyperarc's final local value across ranks that hold pieces of $h$, yielding the global subtree measure
$$
S(h)=\sum_r s_r(h,\mathrm{end}).
$$

For a hyperarc with ordered sequence $a_1 \to a_2 \to \cdots \to a_k$, the local subtree contribution on rank $r$ is
$$
s_r(i)=\sum_{j=1}^i m_r(a_j), \qquad S_r=s_r(k).
$$
The global quantity is then computed as
$$
S(h)=\sum_{r \in \mathrm{owners}(h)} S_r
$$
via `MPI_Allreduce (sum)`. For boundary arcs, designated as ghost arcs, each mesh vertex is assigned exclusively to the rank of highest block-ID that holds it, specifically to avoid double-counting.

The complexity summary is explicit. Local prefix scans cost $O(|A_r|/\mathrm{threads} + \log |h|)$ time. `MPI_Allreduce` on $O(|H_r|)$ hyperarcs costs $O(\log P)$ latency plus $O(|H_r|)$ words. Overall time per hypersweep is approximately
$$
O((t/P) + \log P).
$$

This formulation situates the DHCT not merely as a topological skeleton but as a distributed query structure for subtree volumes and related geometric properties. A plausible implication is that the ordered hyperarc representation is the critical bridge between topological connectivity and scalable prefix-style aggregation.

## 4. Branch decomposition and significance measures

The branch decomposition is computed in three stages: assigning each supernode its best ascending and descending continuation, exchanging those choices to ensure global consistency, and applying pointer-doubling to identify each monotone branch set [2408.04836].

In the first stage, each supernode $v \in V_r$ examines its up-adjacent arcs $\{u \to v\}$ and down-adjacent arcs $\{v \to w\}$ and selects the one whose significance $\phi(a)$ is maximal. The report states that significance may be arc-length, estimated volume, or persistence. In the second stage, over the same $R$ rounds used for augmentation fan-in, partner ranks exchange, for each shared supernode, $\mathrm{best\_up}(v)$ and $\mathrm{best\_down}(v)$. Upon receipt, each rank updates
$$
\mathrm{best\_up}(v) := \arg\max \{ \phi(\mathrm{best\_up}(v)), \phi(\mathrm{best\_up}'(v)) \},
$$
and similarly for $\mathrm{best\_down}$. After $R$ rounds, all ranks agree on these choices for every shared supernode.

In the third stage, each chosen best-up pointer is oriented toward the root of the hyperstructure, and pointer-doubling, described as parallel tree-contraction, is used to find the root representative $r_v$ for each vertex in $O(\log \mathrm{depth})$ steps. The set
$$
\{ v \mid r_v = c \}
$$
forms a branch $B_c$ rooted at $c$.

The significance metrics are defined directly from arc and hyperarc quantities. For an arc $a=(u \to v)$ with end isovalues $f(u)>f(v)$,
$$
\mathrm{height}(a)=|f(u)-f(v)|, \qquad \mathrm{volume}(a)=S(h_a),
$$
where $S(h_a)$ is the subtree volume obtained from the hypersweep on hyperarc $h_a$. The persistence of a branch $B$ is
$$
\mathrm{pers}(B)=\max_{a \in B} \mathrm{height}(a),
$$
and its volume can be written as
$$
\mathrm{vol}(B)=\sum_{a \in B} \mathrm{volume}(a).
$$
The report states that $\mathrm{vol}(B)$ is typically used as $\phi(B)$. After pointer-doubling, a local reduction by root-ID computes branch volumes or persistences through `MPI_Allreduce by key (r_v) sum (\phi(a))`.

This decomposition supplies the ranking mechanism needed for later contour queries. It also clarifies that branch extraction is not presented as an independent geometric post-process; instead, it is integrated with the distributed topological structure through shared continuation choices and global reductions.

## 5. Query structure and contour extraction

After branch significances $\phi(B)$ have been computed, the DHCT is used as a query structure for selecting important contours. The selection step chooses the top-$k$ branches $\{B_1,\ldots,B_k\}$, and one representative contour per branch is extracted by choosing an isovalue at the branch's saddle-end:
$$
\eta_B = f(v_s) + \epsilon,
$$
where $v_s$ is the lower, saddle end of $B$, and $\epsilon$ is an infinitesimal value above $f(v_s)$. In practice, the implementation uses the next-higher floating-point value [2408.04836].

The distributed query proceeds in four phases. First, top-$k$ branches are selected globally by repeated partner exchanges over levels $\ell=1,\ldots,R$, merging local top-$k$ lists. Second, for each selected branch, the corresponding saddle supernode $v_s$ is found and $\eta[B]$ is set to `nextfloat(f(v_s))`. Third, each rank extracts local contours by iterating over each cell $c$ in $M_r$ and each selected branch $B$, testing whether $\min(f(c)) < \eta[B] < \max(f(c))$. When the test passes, `MarchingCells` computes a fragment $F(c,B)$. The local hyperarc $h_c$ for that fragment is identified, and its branch identifier $B_h$ is obtained through hyperstructure lookup; the fragment is appended to `local_mesh[B]` only if $B_h = B`. Fourth, the per-branch meshes may optionally be gathered or rendered, with `MPI_Gatherv` used if global mesh assembly is desired.

The communication and local-work characteristics are narrowly specified. In step 1, communication is $O(r\,k)$ words. Step 3 requires no communication because cells remain local. Final gathering is optional. Local work consists of $O(N_r \cdot k)$ cell tests plus $O(\log t)$ hyperstructure lookups each.

This query model emphasizes that the branch decomposition is operational rather than purely descriptive. A common misunderstanding would be to treat topological branches as only abstract summaries; here, they are used to drive explicit contour extraction in distributed memory.

## 6. Performance characteristics and scaling behavior

The reported implementation was evaluated on NERSC Perlmutter using AMD EPYC 7763, $2 \times 64$ cores per node, InfiniBand, up to 512 nodes, 32 768 MPI ranks, and 128 OpenMP threads per node [2408.04836]. The four datasets tested were GTOPO, WarpX, Nyx, and MICrONS.

| Dataset | Single-node time (s) | Speed-up |
|---|---:|---:|
| GTOPO | 85.4 | 17.2× |
| WarpX | 102.1 | 98.8× |
| Nyx | 519.3 | 7.5× |
| MICrONS | — | 6.8× |

The dataset descriptions and headline scaling results are given explicitly. GTOPO is a 2D dataset of size $21\,601 \times 43\,201$, with speed-up approximately $17\times$ at 512 nodes over single-node PPP. WarpX is a 3D dataset of size $6791 \times 371 \times 371$, with speed-up approximately $99\times$ at 256 nodes. Nyx is a 3D dataset of size $1024^3$, with speed-up approximately $7.5\times$ at 512 nodes. MICrONS is 3D $1024^3$ EM data and exhibits scaling similar to Nyx.

The phase-by-phase analysis distinguishes between topological regimes. In strong scaling at fixed problem size, phase (1), DHCT fan-in/fan-out, is dominated by local SMP merge-tree cost and scales well. Phases (4), augmentation, (5), hypersweep, and (6), branch decomposition, drop roughly by one half when doubling nodes for clean topology datasets such as WarpX and GTOPO. For noisy topology datasets such as Nyx and MICrONS, communication of boundary and hyperarc data becomes the bottleneck, flattening strong scaling beyond approximately 128 nodes. Phase (7), contour extraction, scales almost perfectly because only $O(r\,k)=O(\log P)$ exchange is required for small $k$, with $k=10$ given as an example.

The hardware configuration is specified as Perlmutter CPU-only nodes with 64 cores per socket, two sockets, 512 GB RAM, EDR InfiniBand, Cray MPICH for MPI, and OpenMP threading. The summary conclusion is that the distributed DHCT algorithms for augmentation, hypersweeps, branch decomposition, and contour extraction achieve near-ideal scaling on data with well-behaved topology and remain feasible, though communication-bound, for highly complex datasets, enabling interactive-scale contour queries on meshes far larger than any single node's memory [2408.04836].

A plausible implication is that the primary scalability constraint is not the contour-tree abstraction itself but the interaction between topological complexity and boundary-induced communication. That interpretation is consistent with the observed divergence between clean-topology and noisy-topology datasets.

Source: https://www.emergentmind.com/topics/distributed-hierarchical-contour-trees