PYTHON · STEP-BY-STEP LESSON

Connected regions, periodic loops, holes and cavities

I write a small part, check it, and then add the next part.

The calculation, one step at a time

I show the purpose of each notebook cell beside its code and recorded result. The controls move through those results; they do not execute Python in the browser. The notebook remains available for running or changing the calculation.

Loading the lesson…

I construct small masks with known topology before analyzing a simulated microstructure. A mask is a Boolean array: True means occupied. In 2D I use four edge-sharing neighbors. In 3D I use six face-sharing neighbors. Distances are measured in cells. A graph loop, a hole in a filled region and a closed boundary contour are different objects.

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.

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 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: Input values

I retain the default values from the original calculation. The output name locates the files written by this notebook.

Step 3: 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 collections import deque
import numpy as np
from skimage.measure import label, euler_number, find_contours

Step 4: Seam joining and periodic winding

I keep a wrapped pixel address and a lifted, unwrapped position for each visited pixel. A loop displacement divided by the grid shape is a winding vector. Its independent directions determine winding rank. Four edge-sharing neighbors define this 2D graph. A seam crossing alone does not establish wrapping.

Read this complete operation (38 lines)

def periodic_components(mask):
    """Four-neighbour components and independent periodic winding vectors.

    The graph measures connectivity. Its many local cycles must NOT be
    interpreted as holes of the filled pixel region.
    """
    mask = np.asarray(mask, dtype=bool)
    shape = np.array(mask.shape)
    assert mask.ndim == 2 and min(shape) >= 3
    seen = np.zeros(mask.shape, dtype=bool)
    lifts = np.zeros((*mask.shape, 2), dtype=np.int64)
    result = []
    for start in map(tuple, np.argwhere(mask)):
        if seen[start]:
            continue
        seen[start] = True
        queue = deque([start])
        count, basis = (0, [])
        while queue:
            p = queue.popleft()
            count += 1
            for step in ((1, 0), (-1, 0), (0, 1), (0, -1)):
                q = tuple((np.array(p) + step) % shape)
                if not mask[q]:
                    continue
                proposed = lifts[p] + step
                if not seen[q]:
                    seen[q] = True
                    lifts[q] = proposed
                    queue.append(q)
                else:
                    difference = proposed - lifts[q]
                    assert np.all(difference % shape == 0)
                    w = difference // shape
                    if np.any(w) and (not basis or (len(basis) == 1 and basis[0][0] * w[1] - basis[0][1] * w[0] != 0)):
                        basis.append(tuple(map(int, w)))
        result.append({'pixels': count, 'winding_rank': len(basis), 'basis': basis})
    return result

Step 5: My check of this result

I check finite values, conservation or the known reference result as specified below. An assertion stops the notebook if the condition fails. Inspect the quantities and tolerance rather than suppressing a failure.

seam = np.zeros((8, 8), bool)
seam[[0, 7], 3:5] = True
assert label(seam, connectivity=1).max() == 2
assert periodic_components(seam)[0]['winding_rank'] == 0
stripe = np.zeros((8, 8), bool)
stripe[:, 3:5] = True
assert periodic_components(stripe)[0]['winding_rank'] == 1
staircase = np.zeros((8, 8), bool)

Step 6: Repeat the operation and inspect each case

I repeat the operation over the specified cases. Read the loop variable first, then follow what changes on each pass. The assertions check the stated numerical or topological condition.

for i in range(8):
    staircase[i, i] = staircase[(i + 1) % 8, i] = True

Step 7: My check of this result

assert periodic_components(staircase)[0]['winding_rank'] == 1
full = np.ones((8, 8), bool)
assert periodic_components(full)[0]['winding_rank'] == 2
diagonal = np.eye(3, dtype=bool)
assert label(diagonal, connectivity=1).max() == 3
assert label(diagonal, connectivity=2).max() == 1
disk = np.zeros((9, 9), bool)
disk[2:7, 2:7] = True
ring = disk.copy()

Step 8: My check of this result

ring[3:6, 3:6] = False
assert euler_number(disk, connectivity=1) == 1
assert euler_number(ring, connectivity=1) == 0
contours = find_contours(ring.astype(float), 0.5, fully_connected='low')
assert len(contours) == 2
assert all((np.allclose(path[0], path[-1]) for path in contours))
ball = np.zeros((11, 11, 11), bool)
ball[2:9, 2:9, 2:9] = True
shell = ball.copy()

Step 9: My check of this result

shell[3:8, 3:8, 3:8] = False
tunnel = ball.copy()
tunnel[4:7, 4:7, :] = False
assert euler_number(ball, connectivity=1) == 1
assert euler_number(shell, connectivity=1) == 2
assert euler_number(tunnel, connectivity=1) == 0

Step 10: Repeat the operation and inspect each case

I repeat the operation over the specified cases. Read the loop variable first, then follow what changes on each pass. The assertions check the stated numerical or topological condition.

for name, region in [('seam island', seam), ('stripe', stripe), ('staircase', staircase), ('full periodic grid', full)]:
    print(name, periodic_components(region))

seam island [{'pixels': 4, 'winding_rank': 0, 'basis': []}]
stripe [{'pixels': 16, 'winding_rank': 1, 'basis': [(-1, 0)]}]
staircase [{'pixels': 16, 'winding_rank': 1, 'basis': [(-1, -1)]}]
full periodic grid [{'pixels': 64, 'winding_rank': 2, 'basis': [(-1, 0), (0, -1)]}]

Step 11: 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.

print('Bounded 2D Euler: disk=1, annulus=0; annulus boundary contours=2')
print('Bounded 3D Euler: ball=1, shell=2, solid with one tunnel=0')
print('All worked examples passed. Ordinary skimage calls above are NOT periodic.')

Bounded 2D Euler: disk=1, annulus=0; annulus boundary contours=2
Bounded 3D Euler: ball=1, shell=2, solid with one tunnel=0
All worked examples passed. Ordinary skimage calls above are NOT periodic.

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 →