Comparison¶

Created: 31 October 2025 Last Modified: 25 November 2025

---

Description¶

This notebook executes and compares two 3D image segmentation pipelines:

• A custom-built pipeline
• Standard library-based methods

It evaluates performance and visual results using synthetic 3D data (overlapping spheres).

Requirements¶

• Python 3.11
• Required libraries:
  • numpy
  • scipy
  • scikit-image (skimage)

In [1]:

# Import required libraries
import time

import numpy as np
from scipy import ndimage as ndi
from skimage import filters, measure, morphology
from skimage.segmentation import watershed

from src.image_processing.analytical_information import compute_volume_and_com
from src.image_processing.distance_transform import chamfer_distance_transform
from src.image_processing.local_extrema import find_local_minima

# Functions for image processing
from src.image_processing.otsu_method import otsu_threshold
from src.image_processing.watershed_segmentation import watershed_3d

# Import custom modules
# Functions to create synthetic 3D data
from src.shape_creation import create_n_spheres_example

# Functions for visualisations
from src.visualisation import (
    plot_2d_slice_with_values,
    plot_3d_volume_voxels,
    plot_panels,
)

# Integrate figures into notebook
%matplotlib notebook
%matplotlib inline
import matplotlib

# Integrate figures into notebook
matplotlib.pyplot.close()
%matplotlib inline

# Import required libraries import time import numpy as np from scipy import ndimage as ndi from skimage import filters, measure, morphology from skimage.segmentation import watershed from src.image_processing.analytical_information import compute_volume_and_com from src.image_processing.distance_transform import chamfer_distance_transform from src.image_processing.local_extrema import find_local_minima # Functions for image processing from src.image_processing.otsu_method import otsu_threshold from src.image_processing.watershed_segmentation import watershed_3d # Import custom modules # Functions to create synthetic 3D data from src.shape_creation import create_n_spheres_example # Functions for visualisations from src.visualisation import ( plot_2d_slice_with_values, plot_3d_volume_voxels, plot_panels, ) # Integrate figures into notebook %matplotlib notebook %matplotlib inline import matplotlib # Integrate figures into notebook matplotlib.pyplot.close() %matplotlib inline

Step 1: Create Synthetic 3D Image¶

To evaluate the segmentation pipelines, we use synthetic 3D data consisting of two spheres embedded in a 3D volume. This controlled setup allows us to test how well each method performs under different spatial configurations.

The create_two_spheres_example function generates a binary 3D image with two spheres. You can adjust the center coordinates and radii of the spheres to simulate different scenarios, such as:

• Clearly separated spheres
• Touching spheres
• Overlapping spheres

This flexibility helps assess how robust each segmentation approach is when dealing with varying degrees of object separation.

In [2]:

# Generate synthetic 3D image with two spheres
# Default settings create two spheres that touch at one surface
# You can adjust the centers and radii to test different scenarios

# Two overlapping spheres example
# image3d = create_n_spheres_example(
#     centres=[(6, 4, 5),(4, 2, 3)], radii=[2,2], intensities=[180,200],volume_size=10
# )

#image3d = create_n_spheres_example(
#    centres=[(8, 9, 8),(7, 6, 4)], radii=[2,3], intensities=[180,200],volume_size=15)

image3d = create_n_spheres_example(
    centres=[(8, 9, 8),(7, 5, 4),(3,3,3)], radii=[3,4,2], intensities=[180,190,140],volume_size=16)

# Generate synthetic 3D image with two spheres # Default settings create two spheres that touch at one surface # You can adjust the centers and radii to test different scenarios # Two overlapping spheres example # image3d = create_n_spheres_example( # centres=[(6, 4, 5),(4, 2, 3)], radii=[2,2], intensities=[180,200],volume_size=10 # ) #image3d = create_n_spheres_example( # centres=[(8, 9, 8),(7, 6, 4)], radii=[2,3], intensities=[180,200],volume_size=15) image3d = create_n_spheres_example( centres=[(8, 9, 8),(7, 5, 4),(3,3,3)], radii=[3,4,2], intensities=[180,190,140],volume_size=16)

Visual Inspection of the Synthetic Data¶

To verify that the synthetic 3D image was generated correctly, we visualize a central slice through the volume. This slice is taken from the middle of the Z-dimension, which typically intersects both spheres — especially when they are touching or overlapping.

This quick visual check helps confirm:

• The relative positions of the two spheres
• Whether they are overlapping, touching, or clearly separated
• That the intensity values are as expected for segmentation

We use the plot_panels function with a single panel for consistency with later visual comparisons.

In [3]:

# Visualize a single 2D slice of the synthetic 3D image using plot_panels
plot_panels(
    n=1,
    data_list=[image3d],
    plot_func=plot_2d_slice_with_values,
    plot_kwargs_list=[{"slice_index": image3d.shape[0] // 2}],
    title="Synthetic 3D Image Slice",
    subtitles=["Central Z-slice"]
)

# Visualize a single 2D slice of the synthetic 3D image using plot_panels plot_panels( n=1, data_list=[image3d], plot_func=plot_2d_slice_with_values, plot_kwargs_list=[{"slice_index": image3d.shape[0] // 2}], title="Synthetic 3D Image Slice", subtitles=["Central Z-slice"] )

Step 2: Thresholding the Image¶

To prepare the synthetic 3D image for segmentation, we need to separate foreground objects (the spheres) from the background. This is typically done by applying a threshold to the image intensity values.

Thresholding can be approached in different ways:

• A manual threshold can be selected by inspecting the image and choosing a value that visually separates the objects.
• Alternatively, a systematic method like Otsu's thresholding can be used. Otsu's algorithm automatically determines an optimal threshold by minimizing the variance within each class (foreground and background).

In this step, we apply both a custom implementation of Otsu's method and the standard version from scikit-image to compare their results.

In [4]:

# Apply Otsu thresholding using both custom and library implementations
otsu_value = otsu_threshold(image3d)  # Custom
otsu_value_lib = filters.threshold_otsu(image3d)  # Library implementation

# Print the computed threshold values for comparison
print(f"Otsu threshold value (Custom): {otsu_value}")
print(f"Otsu threshold value (Library): {otsu_value_lib}")

# Create binary masks using the computed thresholds
# A binary image is a simplified version of the original where each voxel is either:
# - True (1): foreground — part of the object (e.g., sphere)
# - False (0): background — everything else
# This binary representation is essential for downstream processing like
# distance transforms and segmentation.

binary = image3d > otsu_value
binary_lib = image3d > otsu_value_lib

# Apply Otsu thresholding using both custom and library implementations otsu_value = otsu_threshold(image3d) # Custom otsu_value_lib = filters.threshold_otsu(image3d) # Library implementation # Print the computed threshold values for comparison print(f"Otsu threshold value (Custom): {otsu_value}") print(f"Otsu threshold value (Library): {otsu_value_lib}") # Create binary masks using the computed thresholds # A binary image is a simplified version of the original where each voxel is either: # - True (1): foreground — part of the object (e.g., sphere) # - False (0): background — everything else # This binary representation is essential for downstream processing like # distance transforms and segmentation. binary = image3d > otsu_value binary_lib = image3d > otsu_value_lib

Otsu threshold value (Custom): 0
Otsu threshold value (Library): 0

Visual Comparison of Binary Masks¶

To verify the effectiveness of the thresholding step, we visualize the binary masks produced by both the custom and library implementations of Otsu's method. This helps confirm that the spheres are correctly segmented and highlights any differences in how the two methods interpret the image.

In [5]:

# Visualize binary masks from custom and library thresholding
slice_index = 7
plot_panels(
    n=2,
    data_list=[binary, binary_lib],
    plot_func=plot_2d_slice_with_values,
    plot_kwargs_list=[
        {"slice_index": slice_index},
        {"slice_index": slice_index}
    ],
    title="Binary Masks from Otsu Thresholding",
    subtitles=["Custom Otsu", "Library Otsu"]
)

# Visualize binary masks from custom and library thresholding slice_index = 7 plot_panels( n=2, data_list=[binary, binary_lib], plot_func=plot_2d_slice_with_values, plot_kwargs_list=[ {"slice_index": slice_index}, {"slice_index": slice_index} ], title="Binary Masks from Otsu Thresholding", subtitles=["Custom Otsu", "Library Otsu"] )

Step 3: Compute Distance Transform¶

Once we have a binary image, we can compute a distance transform, which measures the distance from each foreground voxel to the nearest background voxel. This is a crucial step in many segmentation pipelines, especially when preparing for watershed segmentation.

There are two common approaches:

• Exact Euclidean Distance Transform: This method calculates the true shortest distance from each foreground voxel to the nearest background voxel using Euclidean geometry. It provides precise results and is ideal when accuracy is critical. In our case, we use a library implementation (scipy.ndimage.distance_transform_edt) that leverages highly optimized, compiled C code under the hood — making it both fast and reliable despite its complexity.

• Chamfer Distance Transform: This is a fast, approximate method that uses weighted local operations to estimate distances. It’s simpler to implement and easier to understand, making it a good candidate for custom code. While less precise than Euclidean, it is significantly faster and often sufficient for many segmentation tasks.

In this step, we apply both methods to the binary image and compare their outputs and performance. This helps us understand the trade-offs between computational cost, implementation complexity, and accuracy when choosing a distance transform for 3D image segmentation.

In [6]:

# Compute chamfer distance transform using custom implementation
# We wrap the calls in a timer to measure execution time
start_time = time.time()
distance_chamfer = chamfer_distance_transform(binary)
elapsed_chamfer = time.time() - start_time
print(f"Chamfer distance transform time: {elapsed_chamfer:.4f} seconds")

# Compute exact Euclidean distance transform using scipy (library implementation)
start_time = time.time()
distance_lib = ndi.distance_transform_edt(binary_lib)
elapsed_lib = time.time() - start_time
print(f"Exact Euclidean distance transform time: {elapsed_lib:.4f} seconds")

# Compute chamfer distance transform using custom implementation # We wrap the calls in a timer to measure execution time start_time = time.time() distance_chamfer = chamfer_distance_transform(binary) elapsed_chamfer = time.time() - start_time print(f"Chamfer distance transform time: {elapsed_chamfer:.4f} seconds") # Compute exact Euclidean distance transform using scipy (library implementation) start_time = time.time() distance_lib = ndi.distance_transform_edt(binary_lib) elapsed_lib = time.time() - start_time print(f"Exact Euclidean distance transform time: {elapsed_lib:.4f} seconds")

Chamfer converged after 4 iterations. 
Chamfer distance transform time: 0.0037 seconds
Exact Euclidean distance transform time: 0.0017 seconds

Visual Comparison of Distance Transforms¶

Now we visualize the distance maps produced by both methods. In this synthetic example, the Chamfer approximation closely resembles the Euclidean result. However, subtle differences appear at higher distance values near the center of the spheres, where the chamfer method's discrete steps become more noticeable.

This confirms that while Chamfer is often sufficient and easy to implement, Euclidean distance offers smoother and more accurate gradients — which can be useful for downstream tasks like watershed segmentation.

In [7]:

# Visualize the distance transform results using a central slice
# Choose a slice for visualization
slice_index = 7
plot_panels(
    n=2,
    data_list=[distance_chamfer,distance_lib],
    plot_func=plot_2d_slice_with_values,
    plot_kwargs_list=[
        {"slice_index": slice_index},
        {"slice_index": slice_index}
    ],
    title="Distance Transform Comparison",
    subtitles=["Chamfer Approximation", "Exact Euclidean"]
)

# Visualize the distance transform results using a central slice # Choose a slice for visualization slice_index = 7 plot_panels( n=2, data_list=[distance_chamfer,distance_lib], plot_func=plot_2d_slice_with_values, plot_kwargs_list=[ {"slice_index": slice_index}, {"slice_index": slice_index} ], title="Distance Transform Comparison", subtitles=["Chamfer Approximation", "Exact Euclidean"] )

Performance Note¶

Interestingly, despite computing exact distances, the Euclidean method is faster in this case. This is due to its implementation in optimized, compiled C code via scipy.ndimage. It highlights the power of using well-optimized scientific libraries — even for computationally intensive tasks.

Step 4: Detect Local Minima¶

To prepare for watershed segmentation, we need to identify markers that represent the core of each object. A common strategy is to find local minima in the distance transform — these are the points furthest from the background, typically near the center of each object.

We use two approaches to detect these minima:

• Library-based method: skimage.morphology.local_minima provides a fast and reliable way to detect local minima. It uses efficient internal logic and is well-suited for general-purpose use.

• Custom method: We also implement a function find_local_minima(img) that explicitly checks each voxel against its 26 neighbors in 3D space. This method ensures that only absolute local minima are selected — voxels that are strictly smaller than all their neighbors. It’s more transparent and educational, though less optimized than the library version.

This dual approach allows us to compare the behavior and precision of both methods, and gives us flexibility in choosing the most appropriate marker strategy for watershed segmentation.

In [8]:

# Custom implementation: find absolute local minima using 26-neighbor comparison
local_minima = find_local_minima(-distance_chamfer)

# Library implementation: use skimage's built-in local minima detection
local_minima_lib = morphology.local_minima(-distance_lib)

# Count and compare the number of detected minima
print(
    f"Number of local minima (Custom): {np.shape(np.argwhere(local_minima))[0]}"
    f" | (Library): {np.shape(np.argwhere(local_minima_lib))[0]}"
)

# Custom implementation: find absolute local minima using 26-neighbor comparison local_minima = find_local_minima(-distance_chamfer) # Library implementation: use skimage's built-in local minima detection local_minima_lib = morphology.local_minima(-distance_lib) # Count and compare the number of detected minima print( f"Number of local minima (Custom): {np.shape(np.argwhere(local_minima))[0]}" f" | (Library): {np.shape(np.argwhere(local_minima_lib))[0]}" )

Number of local minima (Custom): 3 | (Library): 3

Step 5: Watershed Segmentation¶

Watershed segmentation is used to separate individual objects by simulating a flooding process from seed points — in this case, the local minima of the distance transform. As the "water" rises, each voxel is assigned to the nearest labeled region, effectively partitioning the image into distinct segments.

We apply watershed using two approaches. One uses a custom-built function that implements a voxel-wise 3D flooding algorithm inspired by Meyer’s method. It processes the image level-by-level, starting from the highest distance values and propagating labels from seed points through the binary mask. Each voxel is compared to its 26 neighbors, and label assignment depends on the configuration of neighboring regions:

• If no neighbors are labeled, the voxel inherits its marker label.
• If all labeled neighbors belong to the same region, it adopts that label.
• If neighbors belong to multiple regions, the voxel is assigned to the region with the most adjacent voxels.

This last case requires hardcoded decision logic to resolve conflicts — for example, choosing the label that appears most frequently among neighboring voxels. While this logic is simple and effective, it highlights the need for manual rule-setting in custom implementations, which can affect consistency and scalability.

The second approach uses skimage.segmentation.watershed, which performs the same task using optimized routines. It accepts the inverted distance map and labeled seed points, and produces a segmented volume where each object is assigned a unique label. This method is fast, robust, and widely used in image processing workflows.

Both implementations rely on the same inputs — the distance map and local minima — and produce comparable results, allowing us to evaluate the trade-offs between precision, control, and performance.

In [9]:

# Custom implementation: Execute the watershed algorithm based on Meyer's method
start_custom = time.time()
labels_custom = watershed_3d(distance_chamfer, binary, local_minima)
end_custom = time.time()

# Extract the final segmentation result from the last flooding level
segmentation_custom = labels_custom[-2]
print(f"Custom watershed completed in {end_custom - start_custom:.2f} seconds.")

# Library implementation: Use skimage's optimized watershed function
start_lib = time.time()
# Label markers
markers_lib = measure.label(local_minima_lib)
segmentation_lib = watershed(-distance_lib, markers_lib, mask=binary_lib)
end_lib = time.time()

print(f"Library watershed completed in {end_lib - start_lib:.2f} seconds.")

# Custom implementation: Execute the watershed algorithm based on Meyer's method start_custom = time.time() labels_custom = watershed_3d(distance_chamfer, binary, local_minima) end_custom = time.time() # Extract the final segmentation result from the last flooding level segmentation_custom = labels_custom[-2] print(f"Custom watershed completed in {end_custom - start_custom:.2f} seconds.") # Library implementation: Use skimage's optimized watershed function start_lib = time.time() # Label markers markers_lib = measure.label(local_minima_lib) segmentation_lib = watershed(-distance_lib, markers_lib, mask=binary_lib) end_lib = time.time() print(f"Library watershed completed in {end_lib - start_lib:.2f} seconds.")

Custom watershed completed in 0.02 seconds.
Library watershed completed in 0.00 seconds.

Visualizing the Segmentation Results¶

To evaluate the outcome of the watershed segmentation, we visualize the labeled regions produced by both methods. Each segmented object is assigned a unique label, which we display using a central slice of the 3D volume.

This comparison allows us to:

• Verify that the segmentation correctly separates individual objects
• Observe how the label boundaries align with the original shapes
• Compare the consistency between the custom and library-based implementations

While both methods aim to achieve the same result, subtle differences may appear due to implementation details — such as how ties between neighboring regions are resolved or how plateaus are handled during flooding.

In [10]:

# Plot 3d voxel render of watershed segmentation
plot_panels(
    n=2,
    data_list=[segmentation_custom, segmentation_lib],
    plot_func=plot_3d_volume_voxels,
    plot_kwargs_list=[
        {"threshold_lo": 1, "threshold_hi": 3},
        {"threshold_lo": 1, "threshold_hi": 3},
    ],
    title=None,  # no overall title in original call
    subtitles=["Build", "Libraries"],
    projection="3d",
)

# Plot 3d voxel render of watershed segmentation plot_panels( n=2, data_list=[segmentation_custom, segmentation_lib], plot_func=plot_3d_volume_voxels, plot_kwargs_list=[ {"threshold_lo": 1, "threshold_hi": 3}, {"threshold_lo": 1, "threshold_hi": 3}, ], title=None, # no overall title in original call subtitles=["Build", "Libraries"], projection="3d", )

Step 6: Analytical Analysis of Segmented Objects¶

This is typically the final step in a scientific image processing pipeline. After preprocessing, segmentation, and labeling, the goal is to extract quantitative information from the image — turning pixels into meaningful measurements.

From the labeled segmentation results, we can compute object-level properties such as:

• Volume: The number of voxels per object, which can be converted to physical units if voxel spacing is known.
• Center of Mass: The spatial centroid of each object, useful for tracking, alignment, or spatial statistics.
• Equivalent Diameter: The diameter of a sphere with the same volume as the object — particularly relevant in synthetic datasets with spherical shapes.

These features are what we ultimately want in most scientific imaging tasks: measurable descriptors that allow us to analyze, compare, and interpret the structures in our data. Whether we're studying cells, particles, or anatomical structures, segmentation is just the means — the real value lies in the measurements we extract from it.

In [11]:

# Extract quantitative information from segmented image
# 1. Volume of each labelled segment in nr. of voxels
# 2. Centre of mass for each labelled region in x, y, and z dimension, in nr. of voxels
labels, volume, CoM = compute_volume_and_com(segmentation_custom)
labels_lib, volume_lib, CoM_lib = compute_volume_and_com(segmentation_lib)

# 3. Equivalent diameter, i.e. the diameter a labelled region would have
# if perfectly spherical
equiv_diameter = ((6 * volume / np.pi) ** (1 / 3))
equiv_diameter_lib = ((6 * volume_lib / np.pi) ** (1 / 3))
print("Custom implementation")
print("Label | Volume (voxels) | COM (z, y, x)           | Eq. Diameter (voxels)")
print("----------------------------------------------------------------------")
for i, lab in enumerate(labels):
    x, y, z = CoM[i]
    com_str = f"({x:.0f}, {y:.0f}, {z:.0f})"  # integers as you requested
    print(f"{lab:5d} | {volume[i]:15d} | {com_str:25s} | {equiv_diameter[i]:.2f}")

print("Library implementation")
print("Label | Volume (voxels) | COM (z, y, x)           | Eq. Diameter (voxels)")
print("----------------------------------------------------------------------")
for i, lab in enumerate(labels_lib):
    x, y, z = CoM_lib[i]
    com_str = f"({x:.0f}, {y:.0f}, {z:.0f})"  # integers as you requested
    print(f"{lab:5d} | {volume_lib[i]:15d} | {com_str:25s} | {equiv_diameter_lib[i]:.2f}")

# Extract quantitative information from segmented image # 1. Volume of each labelled segment in nr. of voxels # 2. Centre of mass for each labelled region in x, y, and z dimension, in nr. of voxels labels, volume, CoM = compute_volume_and_com(segmentation_custom) labels_lib, volume_lib, CoM_lib = compute_volume_and_com(segmentation_lib) # 3. Equivalent diameter, i.e. the diameter a labelled region would have # if perfectly spherical equiv_diameter = ((6 * volume / np.pi) ** (1 / 3)) equiv_diameter_lib = ((6 * volume_lib / np.pi) ** (1 / 3)) print("Custom implementation") print("Label | Volume (voxels) | COM (z, y, x) | Eq. Diameter (voxels)") print("----------------------------------------------------------------------") for i, lab in enumerate(labels): x, y, z = CoM[i] com_str = f"({x:.0f}, {y:.0f}, {z:.0f})" # integers as you requested print(f"{lab:5d} | {volume[i]:15d} | {com_str:25s} | {equiv_diameter[i]:.2f}") print("Library implementation") print("Label | Volume (voxels) | COM (z, y, x) | Eq. Diameter (voxels)") print("----------------------------------------------------------------------") for i, lab in enumerate(labels_lib): x, y, z = CoM_lib[i] com_str = f"({x:.0f}, {y:.0f}, {z:.0f})" # integers as you requested print(f"{lab:5d} | {volume_lib[i]:15d} | {com_str:25s} | {equiv_diameter_lib[i]:.2f}")

Custom implementation
Label | Volume (voxels) | COM (z, y, x)           | Eq. Diameter (voxels)
----------------------------------------------------------------------
    1 |              27 | (3, 3, 3)                 | 3.72
    2 |             252 | (7, 5, 4)                 | 7.84
    3 |             118 | (8, 9, 8)                 | 6.09
Library implementation
Label | Volume (voxels) | COM (z, y, x)           | Eq. Diameter (voxels)
----------------------------------------------------------------------
    1 |              29 | (3, 3, 3)                 | 3.81
    2 |             249 | (7, 5, 4)                 | 7.81
    3 |             120 | (8, 9, 8)                 | 6.12

In [ ]: