Automated Cell Segmentation with scikit-image

Courses

Introduction to Calcium Imaging Analysis

ROI Selection

Authors

Dr. Sangeetha Nandakumar | Dr. Nicholas Del Grosso

Open in JupyterLite Download Exercises

Setup

Import Libraries

import numpy as np
import tifffile
import matplotlib.pyplot as plt
from skimage.filters import gaussian
from skimage.feature import peak_local_max
from skimage.morphology import remove_small_objects, binary_closing, disk
from skimage.measure import label, regionprops
from skimage.color import label2rgb
from skimage.exposure import equalize_adapthist
from skimage.segmentation import watershed
from scipy import ndimage as ndi

Download Data

import owncloud
import os

if not os.path.exists('data'):
    print('Creating directory for data')
    os.mkdir('data')

if not os.path.exists('data/demoMovie.tif'):
    oc = owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/NwtdrIg5wGdeGcB')
    oc.get_file('/', 'data/demoMovie.tif');

Creating directory for data

data = tifffile.imread('data/demoMovie.tif')
movie = (data - data.min(axis=(1, 2), keepdims=True)) / (np.ptp(data, axis=(1,2), keepdims=True))
movie.shape

(2000, 60, 80)

This module introduces a structured approach for identifying neurons automatically in calcium imaging datasets. Manual annotation is useful for small-scale analysis, but for large experiments, automated Region of Interest (ROI) detection is essential. The process will be divided into four stages — enhancing the image, intuition of candidate regions, refining masks, and evaluating results.

Section 1: Section 1: Preprocessing

Raw calcium imaging data often contains noise that can obscure the structure of neurons. Before beginning the detection process, it is important to enhance the image so that relevant features are easier to identify. In this section, you will generate projection images from the full movie and apply Gaussian filtering to reduce high-frequency noise.

Code	Description
`np.mean(movie, axis=0)`	Compute the mean of `movie` across all frames (`axis=0` corresponds to frames).
`gaussian(proj, sigma=1)`	Apply a Gaussian filter to `proj` with a sigma value of 1 to smooth the data.

Exercises

Example: Use gaussian smoothing with sigma=1 on mean projection.

proj = np.mean(movie, axis=0)
proj_smooth = gaussian(proj, sigma=1)
plt.imshow(proj_smooth, cmap='gray')

Exercise: Use gaussian smoothing with sigma=10 on mean projection.

Solution

proj = np.mean(movie, axis=0)
proj_smooth = gaussian(proj, sigma=10.0)
plt.imshow(proj_smooth, cmap='gray')

Exercise: Use gaussian smoothing with sigma=0.5 on mean projection.

Solution

proj = np.mean(movie, axis=0)
proj_smooth = gaussian(proj, sigma=0.5)
plt.imshow(proj_smooth, cmap='gray')

Example: Compare mean projection with smoothed (sigma=0.5)

proj = np.mean(movie, axis=0)
proj_smooth = gaussian(proj, sigma=0.5)

plt.subplot(121)
plt.imshow(proj, cmap='gray')

plt.subplot(122)
plt.imshow(proj_smooth, cmap='gray')

Exercise: Compare max projection with smoothed (sigma=0.5)

Solution

proj = np.max(movie, axis=0)
proj_smooth = gaussian(proj, sigma=0.5)

plt.subplot(121)
plt.imshow(proj, cmap='gray')

plt.subplot(122)
plt.imshow(proj_smooth, cmap='gray')

Exercise: Compare standard deviation projection with smoothed (sigma=0.5)

Solution

proj = np.std(movie, axis=0)
proj_smooth = gaussian(proj, sigma=0.5)

plt.subplot(121)
plt.imshow(proj, cmap='gray')

plt.subplot(122)
plt.imshow(proj_smooth, cmap='gray')

Section 2: Section 2: Segmenting Neurons

In this section, you will define ROI boundaries by applying image segmentation techniques. This includes thresholding the smoothed image, applying morphological operations to clean and close the binary mask.

Code	Description
`image > np.percentile(image, 1)`	Create a binary mask where pixel values are greater than the 1st percentile of the image.
`plt.contour(binary_mask, colors='red')`	Plot the contour of the binary mask in red.
`remove_small_objects(binary_mask, min_size=2)`	Remove small objects from the binary mask that are smaller than 2 pixels.
`binary_closing(cleaned_mask, footprint=disk(3))`	Apply binary closing to the `cleaned_mask` with a disk-shaped footprint of size 3 to fill small gaps.

Exercises

proj = np.std(movie, axis=0)
proj_smooth = gaussian(proj, sigma=0.5)

Example: Highlight image areas above 1 percentile of smoothed projection.

image = proj_smooth.copy()
binary_mask = image > np.percentile(image, 1)
plt.imshow(image, cmap='gray')
plt.contour(binary_mask, colors='red')

Exercise: Highlight image areas above 50 percentile of smoothed projection.

Solution

image = proj_smooth.copy()
binary_mask = image > np.percentile(image, 50)
plt.imshow(image, cmap='gray')
plt.contour(binary_mask, colors='red')

Exercise: Highlight image areas above 90 percentile of smoothed projection.

Solution

image = proj_smooth.copy()
binary_mask = image > np.percentile(image, 90)
plt.imshow(image, cmap='gray')
plt.contour(binary_mask, colors='red')

Example: Clean mask by removing any identified objects smaller than 2px.

cleaned_mask = remove_small_objects(binary_mask, min_size=2)
plt.imshow(image, cmap='gray')
plt.contour(cleaned_mask, colors='red')

Exercise: Clean mask by removing any identified objects smaller than 50px.

Solution

cleaned_mask = remove_small_objects(binary_mask, min_size=50)
plt.imshow(image, cmap='gray')
plt.contour(cleaned_mask, colors='red')

Exercise: Clean mask by removing any identified objects smaller than 5px.

Solution

cleaned_mask = remove_small_objects(binary_mask, min_size=5)
plt.imshow(image, cmap='gray')
plt.contour(cleaned_mask, colors='red')

Example: Apply morphological closing using a disk of radius 3.

closed_mask = binary_closing(cleaned_mask, footprint=disk(3))
plt.imshow(image, cmap='gray')
plt.contour(closed_mask, colors='red')

Exercise: Apply morphological closing using a disk of radius 10.

Solution

closed_mask = binary_closing(cleaned_mask, footprint=disk(10))
plt.imshow(image, cmap='gray')
plt.contour(closed_mask, colors='red')

Exercise: Apply morphological closing using a disk of radius 0.6

Solution

closed_mask = binary_closing(cleaned_mask, footprint=disk(0.6))
plt.imshow(image, cmap='gray')
plt.contour(closed_mask, colors='red')

Section 3: Section 3: Identifying ROIs that need Splitting

Some of the ROIs may encompass more than one neuron. As imaging resolution and signal complexity increase, it becomes more common for automated methods to group nearby neurons into a single ROI.

Code	Description
`labeled_rois == 16`	Create a binary mask where the labeled ROIs are equal to `16`.
`proj * roi`	Multiply the projection (`proj`) by the ROI mask (`roi`) to apply the region of interest (ROI) to the image.
`peak_local_max(roi_crop, min_distance=5)`	Detect local peaks in the cropped ROI (`roi_crop`) with a minimum distance of 5 pixels between peaks.
`equalize_adapthist(roi_crop_norm, clip_limit=0.4)`	Apply adaptive histogram equalization to normalize the contrast in `roi_crop_norm` with a clip limit of `0.4`.

Exercises

labeled_rois = label(closed_mask)
np.unique(labeled_rois)

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16],
      dtype=int32)

plt.imshow(image, cmap='gray')
plt.contour(labeled_rois, colors='red')

for roi_num in np.unique(labeled_rois)[1:]:
    temp_roi = labeled_rois == roi_num
    y, x = np.where(temp_roi)
    plt.text(x[0], y[0], str(roi_num))

Example: Show 16th roi.

roi = labeled_rois == 16
plt.imshow(roi)

Exercise: Show 12th roi.

Solution

roi = labeled_rois == 12
plt.imshow(roi)

Exercise: Show 11th roi.

Solution

roi = labeled_rois == 11
plt.imshow(roi)

Example: Detect local peaks in ROI number 16 using a minimum distance of 5

roi = labeled_rois == 16
roi_crop = proj * roi
coords = peak_local_max(roi_crop, min_distance=5)
plt.imshow(roi_crop, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Detect local peaks in ROI number 12 using a minimum distance of 5

Solution

roi = labeled_rois == 12
roi_crop = proj * roi
coords = peak_local_max(roi_crop, min_distance=5)
plt.imshow(roi_crop, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Detect local peaks in ROI number 11 using a minimum distance of 5

Solution

roi = labeled_rois == 11
roi_crop = proj * roi
coords = peak_local_max(roi_crop, min_distance=5)
plt.imshow(roi_crop, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Example: Normalize the cropped ROI and enhance its contrast using adaptive histogram equalization with a clip limit of 0.4.

roi = labeled_rois == 11
roi_crop = proj * roi

roi_crop_norm = roi_crop / roi_crop.max()
enhanced = equalize_adapthist(roi_crop_norm, clip_limit=0.4)
plt.imshow(enhanced, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Normalize the cropped ROI and enhance its contrast using adaptive histogram equalization with a clip limit of 1.0

Solution

roi_crop_norm = roi_crop / roi_crop.max()
enhanced = equalize_adapthist(roi_crop_norm, clip_limit=1.0)
plt.imshow(enhanced, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Normalize the cropped ROI and enhance its contrast using adaptive histogram equalization with a clip limit of 0.1

Solution

roi_crop_norm = roi_crop / roi_crop.max()
enhanced = equalize_adapthist(roi_crop_norm, clip_limit=0.1)
plt.imshow(enhanced, cmap='gray')
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Section 4: Section 4: Watershed Algorithm For Splitting ROIs

The watershed algorithm is a technique that treats the image like a topographic surface, where pixel intensity represents elevation. It “floods” the image from marked points (called markers), and boundaries (edges) are determined by the “watershed lines” where the flooding from different markers meets.

The watershed algorithm uses markers to start the flooding process and a distance map to guide the flow. A distance map is an image where each pixel represents the distance from that pixel to the nearest background pixel

Markers are just the co-ordinates of local intensity peaks which should be converted to labeled markers.

Code	Description
`ndi.distance_transform_edt(enhanced)`	Compute the Euclidean distance transform of the `enhanced` image, where each pixel’s value is the distance to the nearest background pixel.
`watershed(-distance, markers, mask=roi)`	Apply the watershed segmentation algorithm to the negative distance map, using `markers` to define initial regions and `roi` as the mask.

Exercises

Example: Split ROI 11.

roi = labeled_rois == 11
roi_crop = proj * roi
roi_crop_norm = roi_crop / roi_crop.max()
enhanced = equalize_adapthist(roi_crop_norm, clip_limit=0.1)
coords = peak_local_max(enhanced, min_distance=5)
markers = np.zeros_like(enhanced, dtype=np.int32)
for i, (y, x) in enumerate(coords):
    markers[y, x] = i + 1

distance = ndi.distance_transform_edt(enhanced)
labels = watershed(-distance, markers, mask=roi)
plt.imshow(labels, cmap="nipy_spectral")
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Split ROI 13.

Solution

roi = labeled_rois == 13  # Create a binary mask for regions labeled as 13
roi_crop = proj * roi  # Apply the mask to the projection to crop the region of interest
roi_crop_norm = roi_crop / roi_crop.max()  # Normalize the cropped region to [0, 1]

enhanced = equalize_adapthist(roi_crop_norm, clip_limit=0.1)  # Enhance contrast using adaptive histogram equalization
coords = peak_local_max(enhanced, min_distance=5)  # Find local maxima (peaks) in the enhanced image

markers = np.zeros_like(enhanced, dtype=np.int32)  
for i, (y, x) in enumerate(coords):  
    markers[y, x] = i + 1  # Assign a unique label to each peak in the marker array

distance = ndi.distance_transform_edt(enhanced)  # Compute the Euclidean distance transform of the enhanced image
labels = watershed(-distance, markers, mask=roi)  # Apply the watershed algorithm to segment the regions
plt.imshow(labels, cmap="nipy_spectral")  
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')

Exercise: Split ROI 9.

Solution

roi = labeled_rois == 9
roi_crop = proj * roi
roi_crop_norm = roi_crop / roi_crop.max()
enhanced = equalize_adapthist(roi_crop_norm, clip_limit=0.1)
coords = peak_local_max(enhanced, min_distance=5)
markers = np.zeros_like(enhanced, dtype=np.int32)
for i, (y, x) in enumerate(coords):
    markers[y, x] = i + 1

distance = ndi.distance_transform_edt(enhanced)
labels = watershed(-distance, markers, mask=roi)
plt.imshow(labels, cmap="nipy_spectral")
plt.scatter(coords[:, 1], coords[:, 0], s=10, c='r')