Max-Tree Ladder v0.2

[f ≥ λ] is the upper level set: every pixel whose value is at least λ. In the picture: all land above the water.

CC(X) splits a set of pixels into connected components: groups that touch each other. In the picture: the islands. So CC([f ≥ λ]) is "the list of islands at water level λ".

Pxλ is the island that contains pixel x at level λ. Without the λ, Px means "use x's own value as the water level".

In 1D, "connected" just means "next to each other". In a 2D image a pixel also touches the pixels above and below it (4-connectivity), or also the diagonal ones (8-connectivity). The idea stays the same.

Left: each black line is one island, drawn at the level where it is born and spanning the pixels it covers. The vertical line connects an island to the island that swallows it.

Right: the same tree drawn the way the paper draws it (Fig. 2, Fig. 3a): root on top, leaves at the bottom. Leaves (ochre) are peaks. The root is the whole signal.

A new node appears only when the island changes: it is born, it grows, or two islands merge. If the water drops and nothing changes, there is no new node.

Rule 1. Each pixel belongs to one node: the highest island it reaches. Pixel p2 in example A belongs to the node at level 2, not to the root, even though the root also contains it.

Rule 2. Each node picks one of its pixels as its representative. The paper calls it the canonical element. Here we pick the leftmost one, because that is what the 1D algorithm in Rung 4 produces. (The paper's 2D section picks the bottom-right one. Any fixed rule works, it just has to be the same rule everywhere.)

Rule 3. Pointers: a representative points to the representative of its parent node. Any other pixel points to its own node's representative. The root's representative points to null.

A tree is canonical when every pixel points straight at a representative, never at an ordinary pixel. The paper's Fig. 3b is not canonical, Fig. 3c is. "Flattening" or "canonicalization" means fixing the pointers so this is true.

The paper also defines an order p ≺ q ("p is lower in the tree than q"): lower gray value first, and pixel position to break ties. It gives two rules: (a) parent(p) ≼ p, your parent is never higher than you, and (b) parent(p) is always a representative. You will need this order again on the GPU rungs, because it is what lets many threads agree on a winner.

At each step, compare the new pixel p with r, the representative of the island you are standing on:

UP  f(p) > f(r). You climb. Push r on the stack, because its island is not finished. p becomes the new r.

SAME  f(p) = f(r). Same flat island. parent[p] = r.

DOWN  f(p) < f(r). Islands above p are ending. Pop every stack entry that is ≥ f(p), linking each one under the next. Then either p lands on a level you already know (join it), or p is a brand-new lower island (it becomes the parent of r).

The paper's representative rule. In the 1D rungs we picked the leftmost pixel of a node. From here on we use the paper's rule (§2.2): the representative is the pixel that comes last in reading order (the bottom-right-most one). That is why H represents the node {D, E, H}.

This rule comes from the order p ≺ q: "p is lower than q" when f(p) < f(q), or when the values are equal and p has the larger index. Every pixel is either lower or higher than any other pixel, never "tied". Remember this: in Rung 9, many threads will use it to agree on a winner.

make-set(p): p becomes a tree of its own. find-root(n): follow parent arrows from n until you reach a pixel with no parent. union(p, r): hang the tree whose root is r under p. FLATTEN: a final pass, darkest first, that makes every pixel point straight at a representative.

Why this is hard on a GPU. The sort is global: you cannot sort one tile without knowing every other tile. And the joins must happen in sorted order, one after the other. A GPU wants thousands of independent jobs, not one long ordered queue. Rung 8 removes the sort completely.

find-peak-root(x, λ): walk down from x while the next node is still at level λ or above. It stops on the lowest node that is still at or above λ. In the picture: the biggest island containing x that is still above water at level λ.

find-level-root(x) is find-peak-root with λ = f(x): it finds the representative of x's own node.

EXCHANGE(parent(b), a): write a into parent(b) and hand back the value that was there before. That old value is where the zipping continues.

if b ≺ a then swap: when a and b have the same level they are the same node. The order ≺ decides which one stays the representative, so the result is always the same.

Here is the paper's big move (§4.1). Start with every pixel as its own one-pixel tree. Then call CONNECT on every pair of neighbours, in any order. Merging sorted lists gives the same answer whatever order you merge them in, so the final tree is always correct. No sort is needed.

atomicMax(parent(b), a) (Algorithm 6): in one unbreakable step, write a into parent(b) only if a is a better parent (higher in the order ≺) than what is there, and always return the old value. Succeed or fail, the thread gets back "the next node I must merge with" and keeps going. Several threads can make progress in the same round.

atomicCAS, compare-and-swap (Algorithm 7): write only if the cell still holds the value I saw earlier. If someone changed it, my write fails and I retry from the start. Only one thread wins per conflict, so there are more retries (paper Fig. 6a vs 6b).

Where the idea comes from (§3.3, Algorithm 4): GPU connected-component labelling (Komura, 2015) already joins labels with atomicMin, "always point to the smaller label". A fixed rule that every thread obeys means no cycles and no lost joins. The max-tree version uses the order ≺ as its rule.

This is also why the paper calls the method lock-free: no thread ever waits for a lock, and some thread always makes progress.

GPU atomics only compare plain integers. The order ≺ compares a gray level first, then an index. The trick is to store each parent entry as one 32-bit record, highest bits first, so that comparing the integers compares the parents the right way.

Kernel: one program launched on the GPU; the three stages are three kernels. Thread block: a group of threads that share fast shared memory; here one block per tile. SM: the hardware unit that runs blocks. Warp: 32 threads that execute the same instruction together. Global memory: the big, slow memory every block can see.

Stage 1 in shared memory: tile indices fit in 16 bits, so the packed 32-bit record from Rung 9 works with atomicMax. At the end, the tile is copied to global memory and its indices converted to 32 bits. Stage 2 in global memory: one thread per pair of pixels across a tile border. Stage 3: flatten every pixel.

The 1D optimization (§5.1, Fig. 8): each thread first builds the tree of one 16-pixel column with the stack algorithm from Rungs 4–5. Each thread owns its column, so no atomics are needed. Then the columns are merged with CONNECT. This halves the CONNECT calls. The paper measured local-tree time cut by about 50% and the whole algorithm about 30% faster (Table 2). Fewer threads doing more work each was faster, because fewer atomic writes collided. The opposite design, one thread per pixel, was slower.

With diagonals, each border pixel touches three pixels on the other side instead of one, so the border needs three times as many CONNECT calls. Without a fix, the paper measured the Base algorithm about 2.2× slower (+120%) and the 1D-optimized one 2.5× slower (§6.2).

The fix rests on one observation: an edge between u and v is useless if u and v are already joined by another path that never dips lower than both of them. Removing that edge changes no island at any water level, so the tree stays the same.

Propositions 1 and 2 (8-connectivity) keep at most one edge per border position, chosen by comparing each pixel with its neighbour along the border. Proposition 3 (4-connectivity) drops the straight edge at position k when a neighbouring position is a higher bridge. With the rules applied, each thread makes at most one CONNECT, and 8-connectivity costs about the same as 4-connectivity (Fig. 13).

A trade-off the paper points out: with 4-connectivity, removing edges also removes parallelism. In the extreme, one thread merges a whole branch while the others sit idle.

More gray levels means more nodes and longer paths. Each CONNECT walks further, through scattered places in slow global memory. On 16-bit images the paper measured the Base algorithm about 30 times slower than on 8-bit images. Merging in global memory took 90% of the time, and the L2 cache hit rate fell below 40% (§6.3).

The 32 threads of a warp move together. If one thread must walk a long path while its neighbours finish quickly, the whole warp waits for the slowest one (Fig. 14). On some images fewer than 3 of the 32 threads in a warp were active.

Once the tree exists, many filters are just "remove some nodes". Here every node smaller than a minimum area is merged into its parent: small bright specks vanish, big bright shapes stay with their exact outline. The paper's Fig. 1a removes text noise the same way, and astronomers use it to separate stars from galaxies.

The headline: on a single image, memory transfer included, the GPU is at least 5 times faster than the best parallel CPU algorithm on comparable hardware. On a stream of images, the transfer of the next image overlaps the work on the current one, so only kernel time counts: about 10 times faster.

3D images: build 2D slices, then merge them in depth like tiles. Full 3D 26-connectivity would mean a huge number of CONNECT calls in global memory. Tree of shapes: it is built with a max-tree as its last step, so this GPU algorithm could speed it up. Alpha-tree: closely related to a min-tree on the image's edge graph, so the method looks adaptable.