How I develop a calculation
one part at a time.
I build the Hoshen–Kopelman algorithm one part at a time. Each step introduces a small piece of code, explains its purpose and checks an expected result before I continue.
I start with a small grid that I can check by eye. I retain the complete code alongside the individual steps.
Python notebook ↓ · Julia notebook ↓ · All notebooks, scripts and C checkpoints ↓
Install NumPy and Jupyter for the Python notebook. For Julia, use Julia and the IJulia kernel; the algorithm itself needs no Julia packages. For C, each checkpoint is a complete program: compile it, run its assertions, then move to the next file. The download includes commands and the recorded checks.
The tutorial uses four-edge adjacency in 2D. Labels, occupancy values and cell counts are dimensionless. Physical area requires the grid spacing. The topology playground animates labeling and explains wrapping, holes and surfaces.
Why one cluster can initially receive two labels
I begin with five occupied cells. The first row has two separated cells, but the next row joins them. All five cells belong to one cluster. A cluster is a set of occupied cells connected through shared edges; touching at a corner alone does not connect them.
Occupied cells After the first pass Final cluster labels 1 0 1 1 0 2 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
I scan from left to right, then move to the next row. The cell at the upper right receives label 2 because its connection to label 1 has not yet been visited. When the two arms meet, I record that label 2 belongs to the same cluster as label 1. I do not rewrite every earlier cell at that moment.
The parent list records links between labels
parent[k] stores the next label to follow from label k. These are integer label numbers, not memory addresses. For example, the list [0, 1, 1, 2] gives the following sequence:
label 3 label 2 label 1
3 ───────────────► 2 ───────────────► 1
parent[3] = 2 parent[2] = 1 parent[1] = 1
stop here
Label 1 is the root, also called the representative: it points to itself. The function find_label follows these links and returns that root. Label 0 is reserved for background. This three-label sequence illustrates lookup; the five-cell example needs only labels 1 and 2.
To merge two labels, I first find both roots. If they differ, merge_labels makes the larger root point to the smaller root. For the two arms, this gives parent[2] = 1. If both labels already lead to the same root, no new link is needed.
Pass 1: assign labels and record connections
scan_grid visits every cell once. At an occupied cell, I inspect the cells above and to the left because they have already been visited. The four-edge connection rule still includes all four directions: connections to the right and below are considered when those cells are visited.
| Already visited neighbors | My next operation |
|---|---|
| Neither is occupied | I create a new label that points to itself. |
| One is occupied | I follow its label to the root and use that root. |
| Both are occupied | I find both roots. If they differ, I merge them; otherwise they already identify one cluster. |
Pass 2: resolve the provisional labels
finish_labels scans the grid again. For each occupied cell, I find its root and assign a final cluster number. Surviving roots receive consecutive numbers starting at 1. This replaces provisional grid labels; the code does not perform path compression on the parent list. I then count the cells in each final cluster.
Hoshen–Kopelman, one operation at a time
I show the grid, parent links and the code used at the current step together. Cell coordinates are (row, column), starting at 0, as in the Python code. A dot is an occupied cell awaiting a label; 0 is background.
The first three operations belong to pass 1. Pass 2 resolves the grid labels. A new label needs no lookup, and labels with the same root need no new merge.
Grid state
Solid outline: current cell. Dashed outline: an already visited neighbor. Cell colors indicate labels, not composition.
Parent links
Active algorithm block
These excerpts show the caller and helper used at this step. Blue lines mark the current operation. The complete runnable functions appear below.
Compare provisional and final labels
Provisional labels
Final labels
This example uses nonperiodic boundaries. Periodic connections are introduced later in the code lesson. The ring has one cluster and one enclosed hole; these label counts do not measure holes.
The same sequence in code
The functions below follow the operations shown above: look up a root, merge roots, scan the grid, then resolve the labels. I retain the small checks after each function.
Step 1: My small test grid
I begin with a small array that I can inspect by eye. A value of 1 identifies an occupied cell; 0 identifies background. Rows run downward and columns run to the right. These values and the labels are dimensionless. I count connections through shared edges, so corner contact alone does not connect two cells.
Python
import numpy as np
mask = np.array([[1, 0, 1],
[1, 1, 1],
[0, 0, 0]], dtype=bool)
assert mask.shape == (3, 3)
assert mask.sum() == 5
print(mask.astype(int))
The same part in Julia
mask = Bool[1 0 1; 1 1 1; 0 0 0]
@assert size(mask) == (3, 3)
@assert sum(mask) == 5
println(Int.(mask))Step 2: Finding the root of a label
The first row appears to contain two separate clusters. Later cells connect them. I therefore need a way to record that two labels belong to the same cluster. In the parent list, parent[k] gives the next label to follow. A label that points to itself is the representative. Label 0 is reserved for background. I first test lookup without scanning a grid.
Python
def find_label(parent, label):
while parent[label] != label:
label = parent[label]
return label
parent = [0, 1, 1, 2]
assert find_label(parent, 3) == 1
assert find_label(parent, 1) == 1
print("Label 3 belongs to the cluster represented by label 1.")
The same part in Julia
function find_label(parent, label)
while parent[label] != label
label = parent[label]
end
return label
end
parent = [1, 1, 2]
@assert find_label(parent, 3) == 1
println("Label lookup passed.")Step 3: Merging two labels
I find the representatives before merging labels. Otherwise I could link the wrong entries or form a cycle in the parent list. I make the larger representative point to the smaller one. This gives predictable labels in the worked example. I do not need a more elaborate data structure for this small lesson.
Python
def merge_labels(parent, first, second):
first_root = find_label(parent, first)
second_root = find_label(parent, second)
root = min(first_root, second_root)
parent[max(first_root, second_root)] = root
return root
parent = [0, 1, 2, 3]
merge_labels(parent, 2, 3)
merge_labels(parent, 1, 3)
assert [find_label(parent, k) for k in (1, 2, 3)] == [1, 1, 1]
print("Three labels now identify one cluster.")
The same part in Julia
function merge_labels!(parent, first, second)
first_root = find_label(parent, first)
second_root = find_label(parent, second)
root = min(first_root, second_root)
parent[max(first_root, second_root)] = root
return root
end
parent = [1, 2, 3]
merge_labels!(parent, 2, 3)
merge_labels!(parent, 1, 3)
@assert [find_label(parent, k) for k in 1:3] == [1, 1, 1]
println("Label merging passed.")Step 4: The first grid scan
I visit rows from top to bottom, and cells within each row from left to right. Only the upper and left neighbors have already been visited. If both are background, I create a label. If one is occupied, I use its label. If both are occupied, I merge their labels. Each shared-edge connection is considered when its later cell is visited. The returned labels are provisional: some different numbers may already refer to the same cluster.
Python
def scan_grid(mask):
rows, columns = mask.shape
labels = np.zeros(mask.shape, dtype=int)
parent = [0]
for row in range(rows):
for column in range(columns):
if not mask[row, column]:
continue
above = labels[row - 1, column] if row > 0 else 0
left = labels[row, column - 1] if column > 0 else 0
if above == 0 and left == 0:
new_label = len(parent)
parent.append(new_label)
labels[row, column] = new_label
elif above == 0:
labels[row, column] = find_label(parent, left)
elif left == 0:
labels[row, column] = find_label(parent, above)
else:
labels[row, column] = merge_labels(parent, above, left)
return labels, parent
labels, parent = scan_grid(mask)
assert labels[0, 0] != labels[0, 2]
assert find_label(parent, labels[0, 0]) == find_label(parent, labels[0, 2])
print("Provisional labels:")
print(labels)
The same part in Julia
function scan_grid(mask)
rows, columns = size(mask)
labels = zeros(Int, rows, columns)
parent = Int[]
for row in 1:rows
for column in 1:columns
if !mask[row, column]
continue
end
above = row > 1 ? labels[row - 1, column] : 0
left = column > 1 ? labels[row, column - 1] : 0
if above == 0 && left == 0
new_label = length(parent) + 1
push!(parent, new_label)
labels[row, column] = new_label
elseif above == 0
labels[row, column] = find_label(parent, left)
elseif left == 0
labels[row, column] = find_label(parent, above)
else
labels[row, column] = merge_labels!(parent, above, left)
end
end
end
return labels, parent
end
labels, parent = scan_grid(mask)
@assert labels[1, 1] != labels[1, 3]
@assert find_label(parent, labels[1, 1]) == find_label(parent, labels[1, 3])
println("Provisional labels: ", labels)Step 5: The final labels and cluster sizes
I now replace each occupied cell label by its representative. I renumber the surviving representatives from 1 upward. The size of a cluster is the number of cells with its final label. For a uniform two-dimensional grid, its area is this count multiplied by the cell area. A count has no units; an area has the units of the squared grid spacing. The largest-cluster fraction is the largest count divided by the occupied-cell count, not the total grid-cell count.
Python
def finish_labels(labels, parent):
final = np.zeros(labels.shape, dtype=int)
numbers = {}
for row in range(labels.shape[0]):
for column in range(labels.shape[1]):
label = labels[row, column]
if label == 0:
continue
root = find_label(parent, label)
if root not in numbers:
numbers[root] = len(numbers) + 1
final[row, column] = numbers[root]
sizes = [int(np.count_nonzero(final == k))
for k in range(1, len(numbers) + 1)]
return final, sizes
final, sizes = finish_labels(labels, parent)
assert sizes == [5]
assert np.array_equal(final > 0, mask)
print("Final labels:")
print(final)
print("Cluster sizes:", sizes)
The same part in Julia
function finish_labels(labels, parent)
final = zeros(Int, size(labels))
numbers = Dict{Int, Int}()
for row in axes(labels, 1)
for column in axes(labels, 2)
label = labels[row, column]
if label == 0
continue
end
root = find_label(parent, label)
if !haskey(numbers, root)
numbers[root] = length(numbers) + 1
end
final[row, column] = numbers[root]
end
end
sizes = [count(==(k), final) for k in 1:length(numbers)]
return final, sizes
end
final, sizes = finish_labels(labels, parent)
@assert sizes == [5]
@assert (final .> 0) == mask
println("Final labels: ", final, "; sizes: ", sizes)Step 6: Add periodic connections as a separate step
I keep the bounded scan unchanged. For periodic data I then connect occupied cells at matching positions on opposite edges. The top cell in a column is adjacent to the bottom cell in that same column. The left and right cells in a row are adjacent in the same way. I merge these labels before producing the final counts. The two-cell example has two bounded clusters and one periodic cluster. It is a small island crossing the displayed seam; it does not form a loop around the whole domain.
Python
def join_periodic_edges(labels, parent):
rows, columns = labels.shape
for column in range(columns):
top = labels[0, column]
bottom = labels[rows - 1, column]
if top and bottom:
merge_labels(parent, top, bottom)
for row in range(rows):
left = labels[row, 0]
right = labels[row, columns - 1]
if left and right:
merge_labels(parent, left, right)
def hoshen_kopelman(mask, periodic=False):
mask = np.asarray(mask, dtype=bool)
if mask.ndim != 2 or min(mask.shape) == 0:
raise ValueError("Use a nonempty two-dimensional grid.")
labels, parent = scan_grid(mask)
if periodic:
join_periodic_edges(labels, parent)
return finish_labels(labels, parent)
seam = np.zeros((5, 5), dtype=bool)
seam[0, 2] = seam[-1, 2] = True
assert hoshen_kopelman(seam)[1] == [1, 1]
assert hoshen_kopelman(seam, periodic=True)[1] == [2]
print("Seam island: 2 bounded clusters; 1 periodic cluster.")
The same part in Julia
function join_periodic_edges!(labels, parent)
rows, columns = size(labels)
for column in 1:columns
top, bottom = labels[1, column], labels[rows, column]
if top != 0 && bottom != 0
merge_labels!(parent, top, bottom)
end
end
for row in 1:rows
left, right = labels[row, 1], labels[row, columns]
if left != 0 && right != 0
merge_labels!(parent, left, right)
end
end
end
function hoshen_kopelman(mask; periodic=false)
@assert ndims(mask) == 2 && minimum(size(mask)) > 0
labels, parent = scan_grid(Bool.(mask))
if periodic
join_periodic_edges!(labels, parent)
end
return finish_labels(labels, parent)
end
seam = falses(5, 5)
seam[1, 3] = seam[end, 3] = true
@assert hoshen_kopelman(seam)[2] == [1, 1]
@assert hoshen_kopelman(seam; periodic=true)[2] == [2]
println("Periodic seam check passed.")Step 7: My checks on small grids
I test an empty grid, a filled grid, a corner contact and a ring. A ring and a filled square each have one connected component. Cluster labeling cannot distinguish their hole counts. I also test a rectangular grid so that accidentally exchanging the row and column dimensions cannot remain hidden.
Python
assert hoshen_kopelman(np.zeros((3, 5), bool))[1] == []
assert hoshen_kopelman(np.ones((3, 5), bool))[1] == [15]
assert hoshen_kopelman(np.eye(3, dtype=bool))[1] == [1, 1, 1]
ring = np.ones((5, 5), dtype=bool)
ring[1:4, 1:4] = False
assert hoshen_kopelman(ring)[1] == [16]
rectangle = np.zeros((3, 5), dtype=bool)
rectangle[1, 0] = rectangle[1, -1] = True
assert hoshen_kopelman(rectangle)[1] == [1, 1]
assert hoshen_kopelman(rectangle, periodic=True)[1] == [2]
print("All small-grid checks passed.")
The same part in Julia
@assert hoshen_kopelman(falses(3, 5))[2] == Int[]
@assert hoshen_kopelman(trues(3, 5))[2] == [15]
corner = falses(3, 3)
for k in 1:3
corner[k, k] = true
end
@assert hoshen_kopelman(corner)[2] == [1, 1, 1]
ring = trues(5, 5)
ring[2:4, 2:4] .= false
@assert hoshen_kopelman(ring)[2] == [16]
rectangle = falses(3, 5)
rectangle[2, 1] = rectangle[2, end] = true
@assert hoshen_kopelman(rectangle)[2] == [1, 1]
@assert hoshen_kopelman(rectangle; periodic=true)[2] == [2]
println("All small-grid checks passed.")Step 8: Thresholds and the two phases
The composition c and threshold are dimensionless here. I define the high-composition phase using c >= 0.5 and the other phase using c < 0.5. I measure each separately. The occupied fraction uses all grid cells in its denominator. For an absent phase I report a largest-cluster fraction of zero by convention. In a real simulation I would load the saved field instead of this small example. I do not infer winding, holes or numerical convergence from these counts.
Python
composition = np.array([[0.8, 0.2, 0.7],
[0.9, 0.6, 0.8],
[0.1, 0.2, 0.1]])
for name, phase in [("High composition", composition >= 0.5),
("Low composition", composition < 0.5)]:
labels, sizes = hoshen_kopelman(phase, periodic=True)
occupied = int(phase.sum())
largest_fraction = max(sizes) / occupied if occupied else 0.0
assert sum(sizes) == occupied
print(name, "occupied fraction:", occupied / phase.size,
"clusters:", len(sizes), "largest-cluster fraction:", largest_fraction)
The same part in Julia
composition = [0.8 0.2 0.7; 0.9 0.6 0.8; 0.1 0.2 0.1]
for (name, phase) in [("High composition", composition .>= 0.5),
("Low composition", composition .< 0.5)]
phase_labels, phase_sizes = hoshen_kopelman(phase; periodic=true)
occupied = sum(phase)
largest_fraction = occupied > 0 ? maximum(phase_sizes) / occupied : 0.0
@assert sum(phase_sizes) == occupied
println(name, ": occupied fraction = ", occupied / length(phase),
"; clusters = ", length(phase_sizes), "; largest fraction = ", largest_fraction)
endWhat the algorithm establishes
Hoshen–Kopelman assigns the same final label to cells connected by the specified neighbor rule. It gives component counts and sizes. A periodic edge connection does not, by itself, prove a loop around the domain. A ring and a filled square can each have one component but different hole counts. I use separate winding and topology measurements for those questions.
References and source credits
The algorithm is due to Hoshen and Kopelman (1976). Tobin Fricke explains the union–find formulation and provides a C implementation. The code above is a new small teaching implementation. Older C files used in my analysis retain their existing credits, including Prof. T. A. Abinandanan, my PhD adviser (Abi in the original source comments); Fricke’s separate source retains his copyright and GPL notice. I retain the original credits on historical source files.