Marker Generation (Local Minima/Maxima)¶

Created: 25 November 2025 Last Modified: 08 December 2025

Description¶

This notebook finds local minima in a 3D image (typically the distance transform) to serve as markers (seeds) for watershed segmentation. We use 26‑connected neighborhoods and label each absolute local minimum with a unique integer ID.

Inputs:

• img — 3D array (e.g., distance transform; higher inside object, lower near boundaries)

Outputs:

• labels — 3D int array of same shape as img, where each local minimum voxel is assigned a unique positive label (0 elsewhere)

Part of Pipeline¶

• Previous Step: Calculating the Distance Transform
• Next Step: Watershed Segmentation (propagate labels as seeds)

Requirements¶

• Python 3.14
• Libraries:
  • numpy

In [1]:

# Import required libraries
import numpy as np

# Integrate figures into notebook
%matplotlib notebook
%matplotlib inline

# Import required libraries import numpy as np # Integrate figures into notebook %matplotlib notebook %matplotlib inline

Inputs & assumptions¶

We pass in the distance transform (e.g., distance to background). Local minima in the distance map correspond to thin / boundary‑adjacent points; for watershed, seeds are often placed at local maxima of inside-object distance. Here, we show the minima variant (useful in some pipelines or when distance is inverted).

In [2]:

# Load inputs
distance = np.load(file='../src/data/distance_transform.npy')    # Distance transform

# Load inputs distance = np.load(file='../src/data/distance_transform.npy') # Distance transform

Step 1 — Padding and neighbor offsets¶

We pad with +inf so boundary voxels can be compared safely. Neighbors are the 26 positions around a voxel: all combinations of offsets in {−1, 0, 1} except (0,0,0).

In [3]:

padded = np.pad(distance, pad_width=1, mode="constant", constant_values=np.inf)
center = padded[1:-1, 1:-1, 1:-1]

dx, dy, dz = np.meshgrid([-1, 0, 1], [-1, 0, 1], [-1, 0, 1], indexing="ij")
offsets = np.stack([dx.ravel(), dy.ravel(), dz.ravel()], axis=1)
offsets = offsets[np.any(offsets != 0, axis=1)]

neighbors = np.stack(
    [
        padded[
            1 + ox : 1 + ox + distance.shape[0],
            1 + oy : 1 + oy + distance.shape[1],
            1 + oz : 1 + oz + distance.shape[2],
        ]
        for ox, oy, oz in offsets
    ],
    axis=0,
)
print("Neighbors tensor shape:", neighbors.shape)  # (26, X, Y, Z)

padded = np.pad(distance, pad_width=1, mode="constant", constant_values=np.inf) center = padded[1:-1, 1:-1, 1:-1] dx, dy, dz = np.meshgrid([-1, 0, 1], [-1, 0, 1], [-1, 0, 1], indexing="ij") offsets = np.stack([dx.ravel(), dy.ravel(), dz.ravel()], axis=1) offsets = offsets[np.any(offsets != 0, axis=1)] neighbors = np.stack( [ padded[ 1 + ox : 1 + ox + distance.shape[0], 1 + oy : 1 + oy + distance.shape[1], 1 + oz : 1 + oz + distance.shape[2], ] for ox, oy, oz in offsets ], axis=0, ) print("Neighbors tensor shape:", neighbors.shape) # (26, X, Y, Z)

Neighbors tensor shape: (26, 16, 16, 16)

Step 2 — Strict minima/maxima and labeling¶

A voxel is a local minimum/maximum if it’s strictly less/more than all 26 neighbors. We then assign a unique label per minimum/maximum voxel (IDs start at 1).

In [4]:

# Change the inequality sign "<", ">", to obtain minima/maxima
is_max = np.all(center > neighbors, axis=0)

labels = np.zeros_like(distance, dtype=np.int32)
coords = np.argwhere(is_max)
for i, (x, y, z) in enumerate(coords, start=1):
    labels[x, y, z] = i

# Write seedpoints out
np.save(file='../src/data/seedpoints.npy', arr=labels)

print("Local maxima count:", len(coords))

# Change the inequality sign "<", ">", to obtain minima/maxima is_max = np.all(center > neighbors, axis=0) labels = np.zeros_like(distance, dtype=np.int32) coords = np.argwhere(is_max) for i, (x, y, z) in enumerate(coords, start=1): labels[x, y, z] = i # Write seedpoints out np.save(file='../src/data/seedpoints.npy', arr=labels) print("Local maxima count:", len(coords))

Local maxima count: 2

Exercise: Change the definition in the function to detect local minima instead of maxima. Run both versions (find_local_minima and find_local_maxima) on your distance transform.

Questions to explore:

• Does the count of detected extrema correspond to the number of spheres you created earlier?
• If it doesn’t match, why might that be the case? (Hint: Consider overlapping spheres, noise, and plateaus in the distance map.)
• Can visualizing the distance transform slices help explain the discrepancy