DEV Community

kito2718
kito2718

Posted on • Edited on

BCT2: Baseline Source Code Walkthrough

#kaggle #python #bioinformatics #celltracking

Eye Catch Image
Biohub - Cell Tracking During Development

Abstruct

  • Explains the structure and lazy-loading mechanics of OME-Zarr datasets.
  • Step-by-step code walkthrough of the 3D cell detection(Min-Max normalization + blob_dog)and nearest neighbor spatiotemporal tracking(cdist).
  • Demonstrates how to avoid Kaggle \"Submission Error\"s through explicit data casting.

Introduction and Background

In this article, I will explain the inner workings of the baseline pipeline for the Kaggle Biohub Cell Tracking competition.

Although this baseline script doesn't use machine learning(and achieves a modest score of 0.505), it covers essential concepts such as OME-Zarr chunk loading, scientific image normalization, 3D centroid extraction, and nearest neighbor matching.

Content

1. Understanding OME-Zarr Directory Structures

Before diving into the code, we must understand how input data is structured. In OME-Zarr, a .zarr file is actually a directory containing metadata and raw binary arrays.

dataset.zarr/
├── .zgroup            # Defines that this is a hierarchical group structure
├── .zattrs            # Contains physical dimensions and resolution scaling metadata
├── 0/                 # The highest resolution image array (store['0'])
│   ├── .zarray        # Describes array shape (T,C,Z,Y,X) and chunk compression
│   ├── 0.0.0.0.0      # Chunk binary data (T.C.Z.Y.X index)
│   └── ...
Enter fullscreen mode Exit fullscreen mode

The OME-Zarr specification standardizes how datasets are read:

  • .zarray defines chunk sizes(e.g., [1, 1, 30, 256, 256]). This means the large 5D array (T, C, Z, Y, X) is chopped up into small binary files containing exactly 1 frame, 1 channel, 30 slices, and 256x256 pixels.
  • 0.0.0.0.0 corresponds to the first chunk: t=0, c=0, z=0..29, y=0..255, and x=0..255(0-indexed).
  • When we open a store via zarr.open(zarr_path, mode='r'), the library only loads setting files like .zgroup and .zattrs. The actual image data is not loaded into memory until we slice it(e.g., arr[t, 0, :, :, :]). This lazy-loading technique is crucial for managing system memory.

2. Pipeline Architecture

Here is the flowchart and sequence of our cell tracking pipeline:

Pipeline Flowchart

Pipeline Sequence

3. Step-by-Step Code Walkthrough

Step 3.1: Open Zarr Store

We load the Zarr store in read-only mode, fetching only the metadata. We target the highest resolution array named '0'.

# Open the Zarr store
store = zarr.open(zarr_path, mode='r')

# Reference the highest resolution array
arr = store['0']
Enter fullscreen mode Exit fullscreen mode

Step 3.2: Iterate Time Frames

We extract the total frames T and loop through them. We use tqdm to display a progress bar.

num_frames = arr.shape[0]
prev_nodes_info = []

for t in tqdm(range(num_frames), desc="Frames"):
    # Process each frame sequentially
Enter fullscreen mode Exit fullscreen mode

Note: tqdm will display real-time statistics like: Frames: 23%|███ | 23/100 [00:45<02:30, 1.95s/it] showing the current progress and expected remaining time.

Step 3.3: Slice 3D Image

We slice the current 3D image (Z, Y, X) for frame t. This operation triggers the disk read and decompresses the specific binary chunk file.

has_channels = (arr.ndim == 5)
if has_channels:
    img_3d = arr[t, 0, :, :, :]
else:
    img_3d = arr[t, :, :, :]
Enter fullscreen mode Exit fullscreen mode

Step 3.4: Min-Max Normalization & 3D Cell Detection

Microscope images often have uneven illumination. We normalize the values to 0.0..1.0 to maximize contrast, and run blob_dog(Difference of Gaussians)to detect cell nuclei. If no cells are found, we log statistics(Min/Max/Mean)for debugging.

def detect_cells_3d(image_3d, min_sigma=2, max_sigma=5, threshold=0.05):
    img_min = image_3d.min()
    img_max = image_3d.max()
    if img_max > img_min:
        img_norm = (image_3d.astype(np.float32) - img_min) / (img_max - img_min)
    else:
        img_norm = np.zeros_like(image_3d, dtype=np.float32)

    blobs = blob_dog(img_norm, min_sigma=min_sigma, max_sigma=max_sigma, threshold=threshold)
    if len(blobs) > 0:
        return blobs[:, :3] # Return (z, y, x)
    return np.empty((0, 3))
Enter fullscreen mode Exit fullscreen mode

Step 3.5: Physical Scale Transformation & Nearest Neighbor Tracking

Microscope pixels are anisotropic(Z-slice depth is 1.625 µm, while XY pixels are 0.40625 µm). We scale our pixel coordinates to micrometer physical space. Then, we use cdist to compute a distance matrix and execute a greedy nearest-neighbor matching algorithm.

# Convert pixel coordinates to physical space (micrometers)
scale_zyx = np.array([1.625, 0.40625, 0.40625])
coords_prev_physical = coords_prev * scale_zyx
coords_curr_physical = coords_curr * scale_zyx

# Link cells within a threshold of 15.0 µm
links = track_frame_to_frame(coords_prev_physical, coords_curr_physical, max_distance=15.0)
Enter fullscreen mode Exit fullscreen mode

Step 3.6: Combine Nodes & Edges and Type Cast

After processing all frames, we merge detected nodes and tracking edges into a unified DataFrame. To prevent \"Submission Error\"s due to NaNs or float types, we cast all columns explicitly to int64.

# Combine DataFrames
df_sub = pd.concat([df_nodes, df_edges], ignore_index=True)
df_sub.insert(0, 'id', range(len(df_sub)))

# Cast types rigidly
df_sub = df_sub.astype({
    'id': 'int64',
    'node_id': 'int64',
    't': 'int64',
    'z': 'int64',
    'y': 'int64',
    'x': 'int64',
    'source_id': 'int64',
    'target_id': 'int64'
})
Enter fullscreen mode Exit fullscreen mode

Step 3.7: Export CSV 

df_sub.to_csv("submission.csv", index=False)
Enter fullscreen mode Exit fullscreen mode

Complete Source Code

Here is the full source code for the baseline pipeline, achieving a score of 0.505 on Kaggle:

!pip install --no-index --find-links=/kaggle/input/datasets/aaaa1597/zarr-offline-installation-wheels/zarr_wheels zarr

import os
import glob
import zarr
import numpy as np
import pandas as pd
from skimage.feature import blob_dog
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
from tqdm import tqdm

###########################################
# 1. Setup and Data Path Verification
# Set up the input data paths in the Kaggle environment.

CANDIDATES = [
    "/kaggle/input/biohub-cell-tracking-during-development",
    "/kaggle/input/competitions/biohub-cell-tracking-during-development",
]

ROOT = "/kaggle/input/biohub-cell-tracking-during-development"
for p in CANDIDATES:
    if os.path.exists(os.path.join(p, "test")):
        ROOT = p
        break

TEST_DIR = os.path.join(ROOT, "test")
print(f"Using TEST_DIR: {TEST_DIR}")

test_zarr_paths = glob.glob(os.path.join(TEST_DIR, "*.zarr"))
print(f"Found {len(test_zarr_paths)} test datasets.")
for p in test_zarr_paths:
    print(f"  {os.path.basename(p)}")

###########################################
# 2. 3D Blob Detection and Tracking Implementation
# Detect cells in each time frame and track them between adjacent frames.

def detect_cells_3d(image_3d, min_sigma=2, max_sigma=5, threshold=0.05):
    """Detect cell centroids (Z,Y,X) from a 3D image."""
    img_min = image_3d.min()
    img_max = image_3d.max()
    if img_max > img_min:
        img_norm = (image_3d.astype(np.float32) - img_min) / (img_max - img_min)
    else:
        img_norm = np.zeros_like(image_3d, dtype=np.float32)

    blobs = blob_dog(img_norm, min_sigma=min_sigma, max_sigma=max_sigma, threshold=threshold)
    if len(blobs) > 0:
        return blobs[:, :3]
    return np.empty((0, 3))

def track_frame_to_frame(coords_prev, coords_curr, max_distance=15.0):
    """Perform nearest neighbor matching between adjacent frames."""
    if len(coords_prev) == 0 or len(coords_curr) == 0:
        return []

    dists = cdist(coords_prev, coords_curr)

    links = []
    used_curr = set()
    for i in range(len(coords_prev)):
        js = np.argsort(dists[i])
        for j in js:
            if j not in used_curr and dists[i, j] <= max_distance:
                links.append((i, j))
                used_curr.add(j)
                break
    return links

###########################################
# 3. Pipeline Execution
# Loop through all test datasets to collect node and edge information.

nodes = []
edges = []

print(f"Processing {len(test_zarr_paths)} datasets.")
for zarr_path in test_zarr_paths:
    dataset_name = os.path.basename(zarr_path).replace(".zarr", "")
    print(f"Processing dataset: {dataset_name}")

    node_counter = 1
    store = zarr.open(zarr_path, mode='r')
    arr = store['0']

    has_channels = (arr.ndim == 5)
    num_frames = arr.shape[0]

    prev_nodes_info = []

    for t in tqdm(range(num_frames), desc="Frames"):
        if has_channels:
            img_3d = arr[t, 0, :, :, :]
        else:
            img_3d = arr[t, :, :, :]

        coords = detect_cells_3d(img_3d, min_sigma=2, max_sigma=5, threshold=0.05)

        if len(coords) == 0:
            min_val = np.min(img_3d)
            max_val = np.max(img_3d)
            mean_val = np.mean(img_3d)
            print(f"Warning: 0 cells detected in dataset '{dataset_name}' at frame {t}. "
                  f"Image Stats -> Min: {min_val}, Max: {max_val}, Mean: {mean_val:.4f}")

        curr_nodes_info = []
        for idx, (z, y, x) in enumerate(coords):
            node_id = node_counter
            node_counter += 1

            nodes.append({
                "dataset": dataset_name,
                "row_type": "node",
                "node_id": int(node_id),
                "t": int(t),
                "z": int(round(z)),
                "y": int(round(y)),
                "x": int(round(x)),
                "source_id": -1,
                "target_id": -1
            })
            curr_nodes_info.append((idx, node_id, (z, y, x)))

        if len(prev_nodes_info) > 0 and len(curr_nodes_info) > 0:
            coords_prev = np.array([item[2] for item in prev_nodes_info])
            coords_curr = np.array([item[2] for item in curr_nodes_info])

            scale_zyx = np.array([1.625, 0.40625, 0.40625])
            coords_prev_physical = coords_prev * scale_zyx
            coords_curr_physical = coords_curr * scale_zyx

            links = track_frame_to_frame(coords_prev_physical, coords_curr_physical, max_distance=15.0)

            if len(links) == 0:
                print(f"Warning: 0 tracking edges created between frame {t-1} and {t} in dataset '{dataset_name}'. "
                      f"Previous node count: {len(prev_nodes_info)}, Current node count: {len(curr_nodes_info)}.")

            for idx_prev, idx_curr in links:
                src_id = prev_nodes_info[idx_prev][1]
                tgt_id = curr_nodes_info[idx_curr][1]

                edges.append({
                    "dataset": dataset_name,
                    "row_type": "edge",
                    "node_id": -1,
                    "t": -1,
                    "z": -1,
                    "y": -1,
                    "x": -1,
                    "source_id": int(src_id),
                    "target_id": int(tgt_id)
                })

        prev_nodes_info = curr_nodes_info

###########################################
# 4. Submission File Generation and Verification

columns_order = ["dataset", "row_type", "node_id", "t", "z", "y", "x", "source_id", "target_id"]

if len(nodes) == 0:
    df_nodes = pd.DataFrame(columns=columns_order)
else:
    df_nodes = pd.DataFrame(nodes)

if len(edges) == 0:
    df_edges = pd.DataFrame(columns=columns_order)
else:
    df_edges = pd.DataFrame(edges)

df_sub = pd.concat([df_nodes, df_edges], ignore_index=True)
df_sub.insert(0, 'id', range(len(df_sub)))
df_sub = df_sub[["id"] + columns_order]

df_sub = df_sub.astype({
    'id': 'int64',
    'node_id': 'int64',
    't': 'int64',
    'z': 'int64',
    'y': 'int64',
    'x': 'int64',
    'source_id': 'int64',
    'target_id': 'int64'
})

print(f"Total rows: {len(df_sub)}")
df_sub.to_csv("submission.csv", index=False)
print("submission.csv has been successfully generated!")
Enter fullscreen mode Exit fullscreen mode

Japanese Version of This Series

You can read the original Japanese version of this article on Zenn:

Conclusion

Now that we have established a baseline, the next logical step is to improve scores by replacing blob_dog with deep learning models(like Cellpose or StarDist)and optimizing matching algorithms(like the Hungarian method).

I hope this helps!

Top comments (0)