Persistent IndexError when running CSE3D function

[stibossu]

Hello,
I am trying to implement the CSE3D function to my image and segmentation masks.
Unfortunately, I am dealing with a persistent IndexError:

`

IndexError Traceback (most recent call last)
Cell In[52], line 2
1 # This returns a dict containing all the per-metric scores plus the overall quality score.
----> 2 seg_eval_results = CSE3D(img_channels=final_image, metric_mask=combined_masks, PCA_model="3Dv1.6")
4 print("Overall quality score:", seg_eval_results["QualityScore"])
5 for metric, value in seg_eval_results.items(): File ~/.conda/envs/harpy_wsl_cellsegmentationevaluator/lib/python3.10/site-packages/CellSegmentationEvaluator/CSE3D.py:150, in CSE3D(img_channels, metric_mask, PCA_model, img4thresh, voxel_size, disksizes, areasizes)
147 cell_type_labels = cell_type(current_mask, img_channels)
148 else:
149 # get uniformity for nucleus and non-nucleus
--> 150 cell_CV, cell_fraction, cell_silhouette = cell_uniformity(
151 current_mask, img_channels, cell_type_labels 152 )
153 avg_cell_CV = np.average(cell_CV) #was cell_CV[0]
154 avg_cell_fraction = np.average(cell_fraction) #was cell_fraction[0]
File ~/.conda/envs/harpy_wsl_cellsegmentationevaluator/lib/python3.10/site-packages/CellSegmentationEvaluator/functions.py:436, in cell_uniformity(mask, channels, label_list)
434 silhouette.append(silhouette_score(feature_matrix_z, labels))
435 for i in range(c):
--> 436 cluster_feature_matrix = feature_matrix[np.where(labels == i)[0], :]
437 cluster_feature_matrix_z = feature_matrix_z[np.where(labels == i)[0], :]
438 CV_current.append(cell_uniformity_CV(cluster_feature_matrix))
IndexError: index 502 is out of bounds for axis 0 with size 502
`

I don't really understand where it is coming from because I have made the mask IDs in my cell, nuclei and cytoplasm masks contiguous and eliminated mask IDs that are not present in all 3 mask layers. So the zeros in my mask layers are background and the 1,...,N indices are my n mask IDs.
I got the feeling that I am overlooking something small, but I hoping on some advice that would help me further apply your tools.
Best regards,
Stijn Bossuyt

[rfmurphy]

Not sure about the cause. Can you send me one of your image files (and perhaps your python script that is calling CSE) and I will try to reproduce the error?

[stibossu]

I unfortunately can't share any data because of intellectual property reasons.
I could not upload my Jupyter Notebook or python file so I added the code that I used to try perform the cell segmentation evaluation below. The comments in larger font are the most instructive. I hope you are able to find the mistake I made with this code.
Thanks in advance!

`
%load_ext autoreload
%autoreload 2
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all" # multiple outputs per cell (not just last)

Import packages

import sys
import spatialdata as sd
import numpy as np
import cv2 as cv2
import matplotlib.pyplot as plt
import pandas as pd
from CellSegmentationEvaluator.CSE3D import CSE3D
from CellSegmentationEvaluator.single_method_eval_3D import single_method_eval_3D
from skimage.segmentation import relabel_sequential

Input image: sdata is my spatial data object. Clahe is my image layer of interest where I retrieve the raw image array that has dimensions (65, 5000, 5000, dtype=uint16).

cellsegmentationevaluation_path = working_dir + "/CellSegmentationEvaluation/"
Multichannel image is currently (C, Y, X)
final_image = sdata['clahe'].values
Insert a dummy Z dimension of size 1
Image: (C, Y, X) → (C, 1, Y, X)
final_image = final_image[:, np.newaxis, :, :]

Input masks: The raw cell and nuclear masks are (5000, 5000, dtype=int32) arrays that resemble the image dimensions. The arrays contain zeros in non-masked areas and the masked areas are annotated at the location in the array where the cells are with the cell ID of the specific segmented cell. These cell IDs are random ID values but later in the code I will relabel the cell IDs to be contiguous with the "relabel_sequential" function.

Load your raw masks as signed ints
raw_cell_masks = sdata.labels["merged_labels_filtered_v2"].values.astype(np.int32)
raw_nuclear_masks = sdata.labels["nuclear_segmentation"].values.astype(np.int32)

Build raw cytoplasm from the cell masks
raw_cytoplasmic_masks = raw_cell_masks.copy()
raw_cytoplasmic_masks[raw_nuclear_masks > 0] = 0

Only evaluate on “matched” cells. That means zeroing out (backgrounding) any cell in all three channels whose ID doesn’t appear in the nuclear mask. After that, relabel once and you’ll have the exact same set of IDs everywhere.

Filter cis‐nuclear cells only
cell_ids = np.unique(raw_cell_masks)
nuc_ids = np.unique(raw_nuclear_masks)
matched = np.intersect1d(cell_ids, nuc_ids)

mask_keep = np.isin(raw_cell_masks, matched)
zero out everything else in all masks
raw_cell_masks[~mask_keep] = 0
raw_nuclear_masks[~mask_keep] = 0
raw_cytoplasmic_masks[~mask_keep] = 0
Each value in combined_masks[1] (the nucleus channel) must be an integer cell ID also present in combined_masks[0] (the cell channel). Otherwise, label‐mapping inside CSE3D will drop or misassign pixels, leading to mismatched counts.

For CSE3D:

• img_channels should in in the order CZYX
• mask channels assumed to be in the order cell, nucleus, cell-nucleus
• metrics will only be calculated for the first three channels in masks
• img4thresh is either an ndarray containing the image to threshold or the value at which to threshold the summed image
• voxel_size is in cubic microns
• disksizes and areasizes are used for foreground thresholding

Zero-based labeling (0…N–1) was attempted because many NumPy routines are 0-based. However, CSE3D’s internal loops assume positive, contiguous IDs starting at 1, with 0 reserved for background.
calculating how many nuclear-only IDs are present in the nuclear masks
cell_ids = set(np.unique(raw_cell_masks))
nuclei_ids = set(np.unique(raw_nuclear_masks))
extra = nuclei_ids - cell_ids
print("Amount of nuclear-only ID positions:", len(extra))
def filter_nuclear_ids(cell_mask: np.ndarray, nuclei_mask: np.ndarray):
"""
Given two indexed masks of shape (Z,Y,X):
- cell_mask[z,y,x] = integer cell ID (0=background)
- nuclei_mask[z,y,x] = integer nucleus ID (0=background)
Zeroes out any nuclei whose ID doesn’t appear in cell_mask.
Returns a cleaned nuclei_mask (in place).
"""
# Zero out any values ≤ 0
nuclei_mask[nuclei_mask < 1] = 0
# Find the set of valid cell IDs (exclude background=0 if you like)
valid_ids = np.unique(cell_mask)
# Create a boolean mask of where nuclei IDs are NOT (~) in valid_ids
invalid = ~np.isin(nuclei_mask, valid_ids)
# Zero out those invalid nuclear pixels
nuclei_mask[invalid] = 0
return nuclei_mask
In summary, the root cause is an off-by-one mismatch between your mask’s label values and the zero-based row indexing that NumPy requires. By relabeling your mask to zero-based IDs or by iterating only over actual labels, you’ll avoid the “index … out of bounds” error.
CSE3D (and most segmentation‐quality metrics) assume that the mask’s integer labels run contiguously from 0 or 1 up to N, with no gaps. If your input mask has IDs like 5, 12, 37, … there are holes at 6–11, 13–36, etc., and internally it will compute “number of labels” too big and try to index past the end of its feature arrays (hence the out‐of‐bounds error).
Filter out any stray nuclear IDs
clean_nuclear_masks = filter_nuclear_ids(raw_cell_masks, raw_nuclear_masks)
Applying the mask relabeling to all segmentation masks
forward_map is an array of length mask.max()+1 such that
new_id = forward_map[old_id]

inverse_map is an array of length offset + num_cells such that
old_id = inverse_map[new_id]
offset=1
Relabel only the cell mask to zero‐based contiguous IDs
cell_masks_relabeled, cell_masks_fwd_map, cell_masks_inv_map = relabel_sequential(raw_cell_masks, offset=offset)

Use the same forward_map to align nuclei & cytoplasm
nuclear_masks_relabeled = cell_masks_fwd_map[clean_nuclear_masks]

cytoplasm_masks_relabeled = cell_masks_fwd_map[raw_cytoplasmic_masks]

I also tried removing the offset of 1 after the relabel_sequential function but then the background becomes -1 and this also did not work.

Before calling CSE3D, always verify that all three ranges should be identical (e.g. 0 to 572). If one of them still shows 1…573 (or has holes), CSE3D will trip over out-of-bounds indexing.

Even with CellSegmentationEvaluator v1.5.15 installed, the internal cell_uniformity routine still loops over range(c) rather than the actual set of labels, so it sometimes does
feature_matrix[502, …]
even though feature_matrix only has rows 0…501. This happens because sklearn.metrics.silhouette_score insists that cluster labels run contiguously from 0…n_samples − 1, and gaps or labels beyond n_samples−1 trigger an out-of-bounds attempt

Your own relabeling with skimage.segmentation.relabel_sequential(raw, offset=1) is actually correct – it produces {0,1,…,N} with no gaps. It’s only the loop inside CSE3D that isn’t respecting which labels survived filtering.
print("Cell mask IDs:", cell_masks_relabeled.min(), "…", cell_masks_relabeled.max())
print("Nuclear mask IDs:", nuclear_masks_relabeled.min(), "…", nuclear_masks_relabeled.max())
print("Cyto mask IDs:", cytoplasm_masks_relabeled.min(), "…", cytoplasm_masks_relabeled.max())
Stack into CZYX masks (with a dummy Z=1)
combined_masks = np.stack([
cell_masks_relabeled[np.newaxis, :, :],
nuclear_masks_relabeled[np.newaxis, :, :],
cytoplasm_masks_relabeled[np.newaxis, :, :]
], axis=0)

Running the segmentation evaluation

This returns a dict containing all the per-metric scores plus the overall quality score.
seg_eval_results = CSE3D(
img_channels=final_image,
metric_mask=combined_masks,
PCA_model="3Dv1.6"
)

print("Overall quality score:", seg_eval_results["QualityScore"])
for metric, value in seg_eval_results.items():
print(f"{metric:>30} : {value}")
`