Microscopy¶
Optical microscopes use an optical system of lenses and visible light to generate a magnified image of a sample. The maximum attainable magnification in current setups, due to physical limitations, is about 1000x.
While other microscopy techniques with higher magnification exist, our use-cases typically require that the sample, e.g., a bacterial culture, can proliferate while being observed. This is not the case with methods like atomic force microscopy (AFM) and electron microscopy.
Typical image processing problems¶
In the following we assume that an image, or a sequence of images, has been obtained. Processing typically requires to:
- detect the presence of particular objects in an image
- additionally, detect the position if the objects, e.g., via a bounding-box or at the pixel-level (segmentation)
- additionally, track the objects along the sequence of images
In this lecture, we will focus on the segmentation problem.
Idealized Data¶
We will start with idealized data. We assume that objects are disk-shaped, but can potentially overlap.
%%capture
# In case skimage is not already installed, run:
# %pip install scikit-image
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import skimage
from skimage.draw import disk
class Canvas:
x: int
y: int
img: np.ndarray
def __init__(self, x: int, y: int) -> None:
"""Creates a canvas.
x: width
y: length
"""
self.x = x
self.y = y
self.img = np.zeros([x, y], dtype='int64')
def disk(self, x: int, y: int, r: int) -> None:
"""Draws a disk.
x, y: center
r: radius
"""
rr, cc = disk((x, y), r)
self.img[rr, cc] = 1
def show(self) -> None:
plt.imshow(self.img, cmap='gray')
Let's draw some circles. For the sake of challenge, two of them will overlap.
c = Canvas(100, 100)
c.disk(20, 20, 10)
c.disk(35, 35, 15)
c.disk(60, 60, 8)
c.show()
A solution: Distance Maxima and Watershedding¶
We now aim to perform a segmentation of this image. In particular, a so-called segmentation mask has the same dimensions as the original image, with each pixel value being the id of the object it belongs to. Typically, the value 0 is reserved for the background. In our case, 3 more values are to be assigned, corresponding to the 3 circles.
We solve this problem via a commonly used approach:
- Generate a distance field to the background.
- Let the local maxima correspond to the centers of the objects.
- Use the watershedding) algorithm to assign the remaining pixels.
from scipy import ndimage as ndi
# compute distance to 0-valued background
distance = ndi.distance_transform_edt(c.img)
assert isinstance(distance, np.ndarray), "ensure type"
plt.imshow(distance, cmap=plt.cm.hot);
We next use peak_local_max(...)
to find local maxima. The function has a parameter that specifies the locality:
- footprint : ndarray of bools, optional
If provided,
footprint == 1
represents the local region within which to search for peaks at every point inimage
.
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
# show local maxima (https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_watershed.html)
def find_watershed(img: np.ndarray, min_distance: int=10, footprint: np.ndarray=np.ones((10,10))) -> tuple[np.ndarray, int]:
# get distance to background (i.e., value = 0)
distance = ndi.distance_transform_edt(img)
assert isinstance(distance, np.ndarray), "ensure type"
# get local (within footprint) maxima of the distance field
# filter via minimum distance between maxima
coords = peak_local_max(distance, labels=img, min_distance=min_distance, footprint=footprint)
# create segmentation labels of pixels
mask = np.zeros(distance.shape, dtype=bool)
mask[tuple(coords.T)] = True
markers, num_features = ndi.label(mask)
return watershed(-distance, markers, mask=img), num_features
labels, num_features = find_watershed(c.img)
plt.imshow(labels, cmap=plt.cm.nipy_spectral)
print("Number of cells:", num_features)
Number of cells: 3
Noise¶
Let's increase the difficulty by adding some noise.
from numpy import dtype
img = c.img + (1.1 * np.random.random(c.img.shape)).astype("int64")
plt.imshow(img, cmap=plt.cm.gray);
And repeat the above approach.
labels, num_features = find_watershed(img)
plt.imshow(labels, cmap=plt.cm.nipy_spectral)
print("Number of cells:", num_features)
Number of cells: 46
The segmentation and the number of cells is clearly off. Let's do some preprocessing.
from scipy.ndimage import binary_erosion, binary_dilation
# shrink
img_new = binary_erosion(img, iterations=3)
# grow
img_new = binary_dilation(img_new, iterations=3)
plt.imshow(img_new, cmap=plt.cm.nipy_spectral);
labels, num_features = find_watershed(img_new)
plt.imshow(labels, cmap=plt.cm.nipy_spectral)
print("Number of cells:", num_features)
Number of cells: 3
Now the segmentation is correct. Let's move to real data next.
Microscopic data of E. coli¶
We have uploaded an image of E. coli that was obtained with a microfluidic setup, an inverted microscope (40x lense), and a Raspberry Pi + camera. The image data/ecoli.png
is a crop of this image showing 3 bacteria.
from skimage.color import rgb2gray
img = skimage.io.imread("data/ecoli.png")
img = rgb2gray(img)
plt.imshow(img, cmap=plt.cm.gray);
img_bin = np.zeros(img.shape).astype("int64")
mask = (img*255 < 180)
img_bin[mask] = 100
plt.imshow(img_bin, cmap=plt.cm.gray)
img_bin = binary_erosion(img_bin, iterations=1)
img_bin = binary_dilation(img_bin, iterations=1)
img_bin = binary_erosion(img_bin, iterations=1)
img_bin = binary_dilation(img_bin, iterations=1)
plt.imshow(img_bin, cmap=plt.cm.gray)
<matplotlib.image.AxesImage at 0x142685cd0>
labels, num_features = find_watershed(img_bin, min_distance=10)
plt.imshow(labels, cmap=plt.cm.nipy_spectral)
print("Number of cells:", num_features)
Number of cells: 3
Again, a correct segmentation was obtained.
Challenging microscopic data¶
The following uses the dataset by Scherr et al.. Please download it and move it into the main folder.
file_name = "microbeSEG_dataset/test/test/img_000.tif"
img = skimage.io.imread(file_name)
# normalize
img = img / np.max(img) * 255
plt.imshow(img, cmap=plt.cm.gray);
img_bin = np.zeros(img.shape).astype("int64")
mask = (img < 50)
img_bin[mask] = 100
img_bin = binary_erosion(img_bin, iterations=2)
img_bin = binary_dilation(img_bin, iterations=2)
plt.imshow(img_bin, cmap=plt.cm.gray);
On the right center of the image, it is already visible that some of the background has been detected as objects. You can play around with the preprocessing parameters and the threshold of 50
, but no satisfying parameters can be found.
Indeed, subsequent segmentation with watershedding yields wrong results.
labels, num_features = find_watershed(img_bin, min_distance=10, footprint=np.ones([10,10]))
plt.imshow(labels, cmap=plt.cm.nipy_spectral)
print("Number of cells:", num_features)
Number of cells: 352
The number of cells has clearly been overestimated. Playing round with the parameters of find_watershed
may slightly improve the result, but the ground truth is not (easily) obtained.
While including different pre-processing steps, such as normalizing and segmenting on local crops, may solve the problem of a global threshold and yield the correct result, finding the right pipeline can be challenging and time consuming. We explore a different approach in the following section.
Deep learning for segmentation¶
Neural networks (NNs) have shown impressive performance on multiple problems. We will use one to demonstrate the use of NNs for segmenting microscopic images. While several libraries for bacteria / cell segmentation exist, we will use a more general segmentation library in the following. The purpose of this is twofold:
- Currently available specialized libraries are limited in several aspects, making it difficult to adapt them to our needs. Typically one tries them all and picks the best for one's need.
- They come as black boxes, and using them provides little insight into which algorithmic steps are used.
- Several of the libraries are outdated, not running on up-to-date Python libraries.
The classical approach to training NNs to solve segmentation is via supervised learning, training the NN with tuples of an input image img
and a corresponding mask mask
. This requires the availability of labeled masks. Luckily the dataset by Scherr et al. comes with such masks.
While coding your own Pytorch NNs is clearly possible, we skip this step and use the library segmentation_models which comes with a set of standard NN architectures for segmentation as well as the possibility to download pre-trained weights.
For example, a model with the Unet architecture is simply created via:
model = smp.Unet(
encoder_name="resnet34", # choose encoder, e.g. mobilenet_v2 or efficientnet-b7
encoder_weights="imagenet", # use `imagenet` pre-trained weights for encoder initialization
in_channels=1, # model input channels (1 for gray-scale images, 3 for RGB, etc.)
classes=OUT_CLASSES, # model output channels (number of classes in your dataset)
)
import segmentation_models_pytorch as smp
OUT_CLASSES = 1
Here, we will use a different architecture called . We will train our model on the data within microbeSEG_dataset/test/test/
. The following code was adapted from an example segmentation of cars.ipynb).
# adapted from https://github.com/qubvel-org/segmentation_models.pytorch/blob/main/examples/cars%20segmentation%20(camvid).ipynb
import os
import cv2
import re
import torch
import numpy as np
import albumentations as A
import matplotlib.pyplot as plt
from torch.utils.data import DataLoader
from torch.utils.data import Dataset as BaseDataset
from torch.optim import lr_scheduler
import segmentation_models_pytorch as smp
import pytorch_lightning as pl
DATA_DIR = "./microbeSEG_dataset/test/test/"
The files have names composed of a prefix, a _
, and a number with leading 0
s.
Let's code a function that returns the images with a certain prefix. We need this to seperate img
from mask
files; and ensure they are ordered by increasing numbers.
def prefix_files(dir: str, prefix: str) -> list[str]:
pattern = prefix + r'_(\d+)'
img_files = []
for filename in os.listdir(dir):
match = re.match(pattern, filename)
if match:
img_files.append((filename, int(match.group(1))))
# sort by number
img_files.sort(key=lambda x: x[1])
return [os.path.join(dir, f[0]) for f in img_files]
prefix_files(DATA_DIR, "img")
['./microbeSEG_dataset/test/test/img_000.tif', './microbeSEG_dataset/test/test/img_001.tif', './microbeSEG_dataset/test/test/img_002.tif', './microbeSEG_dataset/test/test/img_003.tif', './microbeSEG_dataset/test/test/img_004.tif', './microbeSEG_dataset/test/test/img_005.tif', './microbeSEG_dataset/test/test/img_006.tif', './microbeSEG_dataset/test/test/img_007.tif', './microbeSEG_dataset/test/test/img_008.tif', './microbeSEG_dataset/test/test/img_009.tif', './microbeSEG_dataset/test/test/img_010.tif', './microbeSEG_dataset/test/test/img_011.tif', './microbeSEG_dataset/test/test/img_014.tif', './microbeSEG_dataset/test/test/img_015.tif', './microbeSEG_dataset/test/test/img_016.tif', './microbeSEG_dataset/test/test/img_017.tif', './microbeSEG_dataset/test/test/img_018.tif', './microbeSEG_dataset/test/test/img_019.tif', './microbeSEG_dataset/test/test/img_020.tif', './microbeSEG_dataset/test/test/img_021.tif', './microbeSEG_dataset/test/test/img_022.tif', './microbeSEG_dataset/test/test/img_023.tif', './microbeSEG_dataset/test/test/img_024.tif', './microbeSEG_dataset/test/test/img_025.tif']
class Dataset(BaseDataset):
"""CamVid Dataset. Read images, apply augmentation transformations.
Args:
images_dir (str): path to images folder
class_values (list): values of classes to extract from segmentation mask
augmentation (albumentations.Compose): data transformation pipeline
(e.g. flip, scale, etc.)
"""
CLASSES = [
"cell",
"background",
]
def __init__(
self,
images_dir,
classes=(),
augmentation=None,
dataset_factor: int=10,
):
self.images_fps = prefix_files(images_dir, "img") * dataset_factor
self.masks_fps = prefix_files(images_dir, "mask") * dataset_factor
# convert str names to class values on masks
self.class_values = [self.CLASSES.index(c.lower()) for c in classes]
self.augmentation = augmentation
def __getitem__(self, i):
image = cv2.imread(self.images_fps[i])
mask = cv2.imread(self.masks_fps[i], flags=2)
# extract certain classes from mask (e.g. cells)
masks = [(mask == v) for v in self.class_values]
mask = np.stack(masks, axis=-1).astype("float32")
# invert mask: 1 <-> segment
mask = 1 - mask
# apply augmentation
if self.augmentation is not None:
sample = self.augmentation(image=image, mask=mask)
image, mask = sample["image"], sample["mask"]
return image.transpose(2, 0, 1), mask.transpose(2, 0, 1)
def __len__(self):
return len(self.images_fps)
# helper function for data visualization
def visualize(**images):
"""PLot images in one row."""
n = len(images)
plt.figure(figsize=(16, 5))
for i, (name, image) in enumerate(images.items()):
plt.subplot(1, n, i + 1)
plt.title(" ".join(name.split("_")).title())
if name == "image":
plt.imshow(image.transpose(1, 2, 0))
else:
plt.imshow(image)
plt.show()
Let's inspect one of the image/mask pairs.
dataset = Dataset(DATA_DIR, classes=["cell"])
# check one of them
image, mask = dataset[0]
visualize(
image=image,
cells_mask=mask.squeeze(),
)