Three-dimensional labeling, winding and morphology

The full 256³ calculation is optional. Small exact-mask checks run by default. Set run_full_calculation=True only when enough memory and time are available. I have not rerun the optional full-volume calculation for this notebook conversion.

I measure a saved field without evolving it. c is a dimensionless composition, the phase threshold is 0.5 and dx=0.25 is dimensionless voxel spacing. Six face-sharing neighbors define phase connectivity. Surface area has units of dx² and integrated mean curvature units of dx within this dimensionless model. A simplified display mesh is not used for these measurements.

I develop this calculation in small steps. I run a cell, inspect its result and test it before adding the next part. Run the cells in order in Jupyter. A code cell defines a part of the calculation or carries it out; defining a function does not run it.

I permit free noncommercial teaching, learning and demonstrations with acknowledgment. Research and commercial use require my written permission. Earlier licensed releases retain their permissions. Read the code-use terms.

Copyright (c) 2026 Saswata Bhattacharyya.

Step 1: Prepare the Python environment

I install the packages listed in the next cell once in the same environment as the notebook. In Jupyter or Colab, remove the comment marker from the installation line below if a package is missing. NumPy stores arrays and evaluates numerical operations. Matplotlib displays plots. Pillow reads and writes images. This installation does not change the model parameters.

# %pip install numpy scipy numba scikit-image matplotlib
import numpy as np
assert np.isfinite(np.array([0.0, 1.0])).all()
print("NumPy is ready:", np.__version__)

NumPy is ready: 2.1.3

Step 2: Prepare the next part of the calculation

I define the following operations. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

from pathlib import Path
import sys, json, gzip, time
import numpy as np
from numba import njit
from scipy import ndimage
from skimage.measure import marching_cubes
from skimage.measure._regionprops_utils import EULER_COEFS3D_26

Step 3: Find a union–find root and its displacement

I follow parent links to the representative root. offset stores the displacement along each parent link. I add these displacements, then shorten the parent chain without losing the displacement from each voxel to its root.

@njit(cache=False)
def root_offset(parent, offset, u):
    r = u
    delta = np.zeros(3, np.int32)
    while parent[r] != r:
        delta += offset[r]
        r = parent[r]
    v = u
    prefix = np.zeros(3, np.int32)
    while parent[v] != v:
        nxt = parent[v]
        old = offset[v].copy()
        parent[v] = r
        offset[v] = delta - prefix
        prefix += old
        v = nxt
    return (r, delta)

Step 4: Scan voxels and merge connected labels

I merge occupied predecessor neighbors during the first pass, then join periodic seams. If two voxels already have the same root, a nonzero residual displacement identifies a winding loop. Arrays index voxels, and each index has one parent and a three-component displacement.

@njit(cache=False)
def hk_scan(mask, periodic=True):
    """Scan occupied cells, union predecessor labels, then join periodic seams.
    Offsets encode lifted displacement from a root. Equal-root edges reveal
    winding loops. Background labels are -1. Work arrays use one slot per voxel.
    """
    n = mask.shape[0]
    size = mask.size
    flat = mask.ravel()
    parent = np.arange(size, dtype=np.int32)
    offset = np.zeros((size, 3), np.int32)
    height = np.zeros(size, np.uint8)
    strides = np.array([n * n, n, 1])
    loops = []
    for seam in range(2 if periodic else 1):
        for u in range(size):
            if not flat[u]:
                continue
            coords = np.array([u // (n * n), u // n % n, u % n])
            for axis in range(3):
                if seam == 0:
                    if coords[axis] == 0:
                        continue
                    v = u - strides[axis]
                    step = -1
                else:
                    if coords[axis] != n - 1:
                        continue
                    v = u - (n - 1) * strides[axis]
                    step = 1
                if not flat[v]:
                    continue
                ru, du = root_offset(parent, offset, u)
                rv, dv = root_offset(parent, offset, v)
                difference = du - dv
                difference[axis] += step
                if ru == rv:
                    assert np.all(difference % n == 0)
                    if np.any(difference):
                        loops.append((ru, difference[0] // n, difference[1] // n, difference[2] // n))
                elif height[ru] >= height[rv]:
                    parent[rv] = ru
                    offset[rv] = difference
                    if height[ru] == height[rv]:
                        height[ru] += 1
                else:
                    parent[ru] = rv
                    offset[ru] = -difference
    labels = np.full(size, -1, np.int32)
    for u in range(size):
        if flat[u]:
            labels[u] = root_offset(parent, offset, u)[0]
    records = np.zeros((len(loops), 4), np.int32)
    for i in range(len(loops)):
        records[i, 0] = root_offset(parent, offset, loops[i][0])[0]
        for j in range(3):
            records[i, j + 1] = loops[i][j + 1]
    return (labels.reshape(mask.shape), records)

Step 5: Count components and independent windings

I count voxels attached to each root and divide by the occupied voxel count to obtain the largest fraction of that phase. I retain independent winding vectors for each component. The primary convention is six face-sharing neighbors, including periodic neighbors.

def label_phase(mask, periodic=True):
    if mask.ndim != 3 or len(set(mask.shape)) != 1 or mask.shape[0] < 3:
        raise ValueError('Expected a cubic 3D mask with at least three voxels per axis')
    labels, loops = hk_scan(mask, periodic)
    roots, counts = np.unique(labels[labels >= 0], return_counts=True)
    basis = {int(r): [] for r in roots}
    for record in loops:
        vectors = basis[int(record[0])]
        w = record[1:]
        if len(vectors) < 3 and np.linalg.matrix_rank(np.array(vectors + [w.tolist()], float)) > len(vectors):
            vectors.append(w.tolist())
    total = int(counts.sum())
    order = np.argsort(-counts)
    components = [{'root': int(roots[i]), 'voxels': int(counts[i]), 'winding_rank': len(basis[int(roots[i])]), 'winding_basis': basis[int(roots[i])]} for i in order]
    return (labels, {'components': len(roots), 'phase_voxels': total, 'phase_fraction': total / mask.size, 'largest_fraction_of_phase': int(counts.max()) / total if total else 0, 'clusters': components, 'method': 'Hoshen–Kopelman style scan / union–find; weighted root displacements; per-component winding', 'neighbors': 6, 'periodic': periodic})

Step 6: Measure shortest graph distance from one seed

I visit face-sharing neighbors in breadth-first order. distance counts voxel steps from the seed; it is not physical time. The returned reached count is the number of occupied voxels visited. A disconnected component cannot be reached from this seed.

@njit(cache=False)
def flood_distance(mask, root):
    n = mask.shape[0]
    size = mask.size
    flat = mask.ravel()
    distance = np.full(size, 65535, np.uint16)
    queue = np.empty(size, np.int32)
    head = 0
    tail = 1
    queue[0] = root
    distance[root] = 0
    while head < tail:
        u = queue[head]
        head += 1
        x = u // (n * n)
        y = u // n % n
        z = u % n
        for axis in range(3):
            stride = n * n if axis == 0 else n if axis == 1 else 1
            coordinate = x if axis == 0 else y if axis == 1 else z
            for step in (-1, 1):
                v = u + step * stride
                if coordinate + step < 0:
                    v += n * stride
                if coordinate + step >= n:
                    v -= n * stride
                if flat[v] and distance[v] == 65535:
                    assert distance[u] < 32766
                    distance[v] = distance[u] + 1
                    queue[tail] = v
                    tail += 1
    return (distance.reshape(mask.shape), tail)

Step 7: Calculate a periodic digital Euler characteristic

I classify local voxel configurations and apply scikit-image's Euler lookup coefficients. Periodic convolution includes seam cells. Six-connected foreground is paired with 26-connected background. Digital topology depends on this convention.

def periodic_euler(mask, neighbors=6):
    """Ohser digital Euler LUT as in skimage, with periodic rather than zero edges.
    6-connected foreground uses 26-connected background; 26 reverses this.
    LUT coefficients are from scikit-image (BSD-3-Clause), not historical code.
    """
    config = np.zeros((3, 3, 3), int)
    config[1:, 1:, 1:] = np.array([[[1, 4], [2, 8]], [[16, 64], [32, 128]]])
    code = ndimage.convolve(mask.astype(np.uint8), config, mode='wrap', output=np.uint8)
    hist = np.bincount(code.ravel(), minlength=256)
    coefficients = np.array(EULER_COEFS3D_26)[::-1] if neighbors == 6 else np.array(EULER_COEFS3D_26)
    numerator = int(coefficients @ hist)
    assert numerator % 8 == 0
    return numerator // 8

Step 8: Measure the unsimplified periodic interface

I extract the c=0.5 surface, include seam cells, and weld paired periodic vertices. Surface area is calculated from triangles. I check that each mesh edge belongs to two triangles before interpreting a closed surface. The Euler characteristic is vertices minus edges plus faces. Curvature is a discrete mesh estimate, not a convergence result.

def interface_morphology(c, dx):
    """Full-resolution periodic c=.5 interface. Include seam cells explicitly.
    Mesh area/curvature use actual triangles before any display simplification.
    The oriented normals point out of the high-c phase (checked on test balls).
    """
    n = c.shape[0]
    extended = np.pad(c, ((0, 1),) * 3, mode='wrap')
    vertices, faces, _, _ = marching_cubes(extended, 0.5, allow_degenerate=False, gradient_direction='ascent')
    triangles = vertices[faces].astype(np.float64) * dx
    cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
    lengths = np.linalg.norm(cross, axis=1)
    normals = cross / lengths[:, None]
    area = float(0.5 * lengths.sum())
    coordinates = np.round(np.mod(vertices, n), 6)
    coordinates[coordinates >= n - 1e-06] = 0
    _, inverse = np.unique(coordinates, axis=0, return_inverse=True)
    welded = inverse[faces]
    oriented = np.concatenate([welded[:, [0, 1]], welded[:, [1, 2]], welded[:, [2, 0]]])
    edge_positions = np.concatenate([triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 1], triangles[:, 0] - triangles[:, 2]])
    face_ids = np.tile(np.arange(len(faces)), 3)
    sorted_edges = np.sort(oriented, axis=1)
    order = np.lexsort((sorted_edges[:, 1], sorted_edges[:, 0]))
    keys = sorted_edges[order]
    starts = np.r_[0, 1 + np.flatnonzero(np.any(keys[1:] != keys[:-1], axis=1))]
    counts = np.diff(np.r_[starts, len(keys)])
    closed = bool(np.all(counts == 2))
    answer = {'threshold': 0.5, 'surface_area': area, 'units': 'dimensionless; dx applied', 'vertices': int(len(np.unique(welded))), 'edges': int(len(starts)), 'triangles': int(len(faces)), 'edge_incidence_two': closed, 'boundary_or_nonmanifold_edges': int(np.count_nonzero(counts != 2)), 'periodic_welding_tolerance_voxels': 1e-06, 'display_mesh_used': False}
    if not closed:
        return answer
    a = order[starts]
    b = order[starts + 1]
    coherent = bool(np.all(oriented[a, 0] == oriented[b, 1]) and np.all(oriented[a, 1] == oriented[b, 0]))
    answer['coherent_orientation'] = coherent
    if not coherent:
        return answer
    chi = answer['vertices'] - answer['edges'] + answer['triangles']
    answer['surface_euler'] = chi
    from scipy.sparse import coo_matrix
    from scipy.sparse.csgraph import connected_components
    edges = keys[starts]
    nv = int(inverse.max() + 1)
    graph = coo_matrix((np.ones(len(edges)), (edges[:, 0], edges[:, 1])), shape=(nv, nv))
    components = connected_components(graph, directed=False, return_labels=False)
    answer['surface_components'] = int(components)
    answer['sum_surface_genera'] = int(components - chi // 2)
    assert chi % 2 == 0
    direction = edge_positions[a]
    edge_length = np.linalg.norm(direction, axis=1)
    direction /= edge_length[:, None]
    na = normals[face_ids[a]]
    nb = normals[face_ids[b]]
    theta = np.arctan2(np.einsum('ij,ij->i', direction, np.cross(na, nb)), np.einsum('ij,ij->i', na, nb))
    answer['integrated_mean_curvature_high_c_outward'] = float(0.5 * np.sum(edge_length * theta))
    answer['integrated_gaussian_curvature'] = float(2 * np.pi * chi)
    answer['curvature_convention'] = 'H=(k1+k2)/2; convex high-c ball has positive integrated H; K integral from Gauss–Bonnet'
    return answer

Step 9: Run the complete example

I keep the original algorithm, parameters, seed and output conventions. I inspect the saved diagnostics before interpreting the pictures. A finite run is not evidence of convergence or equilibrium.

def run(field, output, dx=0.25):
    out = Path(output)
    out.mkdir(parents=True, exist_ok=True)
    c = np.load(field)
    n = c.shape[0]
    report = {'field': Path(field).name, 'grid': list(c.shape), 'dx': dx, 'thresholds': {}, 'morphology': {}}
    packed = np.zeros(c.shape, np.uint16)
    flood = {}
    for threshold in (0.45, 0.5, 0.55):
        results = {}
        for phase, mask in [('red', c >= threshold), ('blue', c < threshold)]:
            print('Labeling', threshold, phase, flush=True)
            labels, result = label_phase(mask)
            result['digital_euler_6'] = periodic_euler(mask, 6)
            result['digital_euler_26'] = periodic_euler(mask, 26)
            result['volume'] = int(mask.sum()) * dx ** 3
            results[phase] = result
            if threshold == 0.5:
                assert result['components'] == 1
                center = np.array([n // 2] * 3)
                occupied = np.argwhere(mask)
                root = occupied[np.argmin(np.sum((occupied - center) ** 2, axis=1))]
                del occupied
                distance, reached = flood_distance(mask, int(np.ravel_multi_index(tuple(root), mask.shape)))
                assert reached == result['phase_voxels']
                packed[mask] = distance[mask] + (32768 if phase == 'red' else 0)
                flood[phase] = {'root': root.tolist(), 'reached_voxels': reached, 'maximum_distance': int(distance[mask].max()), 'distance': 'shortest face-sharing steps from one seed, on the full periodic phase graph'}
            del labels
            print(json.dumps(result), flush=True)
        report['thresholds'][str(threshold)] = results
    report['flood'] = flood
    with gzip.open(out / 'phase-flood.bin.gz', 'wb', compresslevel=9) as f:
        f.write(packed.astype('<u2').tobytes())
    print('Measuring full-resolution periodic interface', flush=True)
    report['morphology'] = interface_morphology(c, dx)
    report['morphology']['source'] = 'New implementation inspired by morphological measures studied by Fiałkowski–Aksimentiev–Hołyst; not their recovered C code.'
    (out / 'analysis.json').write_text(json.dumps(report, indent=2))
    print('Morphology', json.dumps(report['morphology']), flush=True)
    return report

Step 10: My check against a known answer

I test the measurement routines on a small exact mask before allocating the full volume. This is a function check, not a smaller replacement for the published 256³ simulation. Notebook-defined Numba functions disable only the disk cache; numerical operations are unchanged.

mask = np.ones((8, 8, 8), dtype=bool)
labels, result = label_phase(mask)
assert result["components"] == 1
assert result["clusters"][0]["winding_rank"] == 3
assert periodic_euler(mask, 6) == 0
shell = np.zeros((8, 8, 8), dtype=bool)
shell[2:6, 2:6, 2:6] = True
shell[3:5, 3:5, 3:5] = False
assert periodic_euler(shell, 6) == 2
surface = interface_morphology(shell.astype(float), 1.0)
assert surface["edge_incidence_two"]
print(result)
print(surface)

{'components': 1, 'phase_voxels': 512, 'phase_fraction': 1.0, 'largest_fraction_of_phase': 1.0, 'clusters': [{'root': 1, 'voxels': 512, 'winding_rank': 3, 'winding_basis': [[0, 0, 1], [0, 1, 0], [1, 0, 0]]}], 'method': 'Hoshen–Kopelman style scan / union–find; weighted root displacements; per-component winding', 'neighbors': 6, 'periodic': True}
{'threshold': 0.5, 'surface_area': 97.40522711209204, 'units': 'dimensionless; dx applied', 'vertices': 120, 'edges': 348, 'triangles': 232, 'edge_incidence_two': True, 'boundary_or_nonmanifold_edges': 0, 'periodic_welding_tolerance_voxels': 1e-06, 'display_mesh_used': False, 'coherent_orientation': True, 'surface_euler': 4, 'surface_components': 2, 'sum_surface_genera': 0, 'integrated_mean_curvature_high_c_outward': 18.84955592153876, 'integrated_gaussian_curvature': 25.132741228718345, 'curvature_convention': 'H=(k1+k2)/2; convex high-c ball has positive integrated H; K integral from Gauss–Bonnet'}

Step 11: The optional full-resolution calculation

I leave this expensive calculation off by default. Setting run_full_calculation to True runs the unchanged original model or full saved-field analysis. Generation and full morphology each need substantial memory. The small tests above do not certify a new full run; the published reports describe the original full-resolution results.

from pathlib import Path
from urllib.request import urlretrieve

work = Path("lesson-data")
work.mkdir(exist_ok=True)

def get_file(relative_path):
    # I download a public example only when its local copy is missing.
    destination = work / Path(relative_path).name
    if not destination.exists():
        urlretrieve("https://saswataiith.github.io/" + relative_path, destination)
    return destination

run_full_calculation = False
if run_full_calculation:
    field_path = get_file("assets/examples/bicontinuity-3d/composition.npy")
    report = run(field_path, "morphology-3d-results", dx=0.25)

Step 12: Review and extend the calculation

I first restart the kernel and run all the required cells. Then I change one input, predict its effect and compare the result with that prediction. A successful check establishes the tested property; it does not establish convergence or equilibrium. I retain the complete script for command-line use. The notebook embeds the calculation directly and does not need that script to run.

Read the notebook execution checks → · More Python, Julia and C lessons →