example cell tracking
This document will walk you through an example of tracking bacteria growing in a microscope time-lapse image.
[[TOC]]
- MATLAB R2016a or later with the Image Processing Toolbox installed
- The Linear Assignment toolbox is installed
- Download the example image file celldemo.tif
celldemo.tif
is a multi-page TIFF file containing time-lapse images of a
growing cyanobacteria cell. We thank the Cameron
Lab for supplying us with the cells to make this
image.
After reviewing this document, you will be able to:
- Write code to segment cells using intensity threshold and watershed algorithms
- Change
LAPLinker
properties to set up tracking - Generate an output image to validate the tracking
- Save the tracked data
- Segmentation - Identification of cells in the image
- Tracking - Identifying and grouping data from the same object between frames
The following are common definitions used to describe object-oriented code:
-
Class, e.g.
LAPLinker
class - The type of object. -
Class file - The m-file that contains the code describing the object, e.g.
LAPLinker.m
-
Object, e.g.
LAPLinker
object - An "instance" or concrete realization of a class. You create an object by assigning a variable name to the class.Example:
obj = LAPLinker
results inobj
which will be aLAPLinker
object. -
Object properties, e.g. LAPLinker properties - Properties contain data for the object. For example, resulting time series data is stored in the property
tracks
of aLAPLinker
object. You access properties of an object using dot notation:obj.tracks
. -
Method - A function defined within the class file that (typically) acts on an object of that class. Methods are called using the object as the first input argument:
treeplot(obj)
.
This section is optional, but recommended if you are new to MATLAB
- Create a new folder in your MATLAB directory and change your working directory
to that folder
mkdir exampleLAPLinker cd exampleLAPLinker
- Download and copy the demo image file into this directory
- Start a new script
edit exampleCells.m
- Copy or type the commands in the following sections into the script
Tip: To get help on a particular function, use
help function_name
. For example,help imread
.
To identify the cells in this image, we will use a simple intensity thresholding algorithm. The basic principle is to identify a threshold intensity that identifies cells in the image (e.g. pixels above or below the threshold are cells and the rest are background).
-
Read and display the first frame of the image using
imwrite
:I = imread('celldemo.tif'); imshow(I, [])
-
Use the datatips tool to determine a suitable threshold.
For this example, since the cells have a dark border, it is easier to set a maximum threshold value. We will choose
1.65e4
for this example. -
Generate a mask (a binary image where true/white pixels indicate a cell and false/black pixels indicate the background)
mask = I < 1.65e4;
-
Display the mask to see the initial results (feel free to run this function after each of the following steps)
imshow(mask, [])
-
Fill in the interior of the cell using
imfill
mask = imfill(mask, 'holes');
-
Use morphological opening to smooth the edges of the cell
mask = imopen(mask, strel('disk', 3));
-
Thicken the cell masks to get a better fit
mask = bwmorph(mask, 'thicken', 3);
-
Remove any objects less than 300 pixels
mask = bwareaopen(mask, 300);
-
Display the final mask to check
imshow(mask)
Once we have a mask, we can use MATLAB's regionprops
function to measure
data about the cell. For this example, we will measure the center coordinate
(centroid
) and the length of the cell (MajorAxisLength
).
For tracking, since these cells are relatively stationary within the image, we
can use the pxintersect
metric to compute the ratio of overlapping pixels to
identify the same object in different frames. To do so, we also need to
get a list of pixels in the cell (PixelIdxList
).
Putting it all together, we get:
data = regionprops(mask, 'Centroid', 'PixelIdxList', 'MajorAxisLength');
As the cells grow and divide, it will be necessary to separate the different cells. For example, try reading in and segmenting frame 10.
I = imread('celldemo.tif', 10);
mask = I < 1.65e4;
mask = imfill(mask, 'holes');
mask = imopen(mask, strel('disk', 3));
imshow(mask)
MATLAB will treat connected regions in the mask as a single object, so we have to split the cells apart. To do so, we will use an algorithm known as the watershed transform. MATLAB has a great explanation of the algorithm here.
We can follow the steps outlined in the documentation for the watershed function:
%Compute the distance transform
dd = -bwdist(~mask);
dd(~mask) = -Inf;
%Perform an h-minima transform to suppress minima with a depth
%less than 1. This will prevent oversegmentation (cell divided
%into too many separate objects)
dd = imhmin(dd, 1);
LL = watershed(dd);
%Update the mask
mask(LL == 0) = 0;
Checking the final mask should show that there are now two objects clearly separated by a line.
imshow(mask)
-
Create a new LAPLinker object
%Create a LAPLinker object L = LAPLinker;
-
Set up tracking parameters
%Set up the tracking parameters to use overlapping pixels L.LinkedBy = 'PixelIdxList'; L.LinkCostMetric = 'pxintersect'; L.LinkScoreRange = [1 12]; L.TrackDivision = true; L.DivisionParameter = 'PixelIdxList'; L.DivisionScoreMetric = 'pxintersect'; L.DivisionScoreRange = [1 12]; L.MinFramesBetweenDiv = 5;
Here we are setting up the object to link objects by the number of overlapping pixels in each frame, as well as tracking division using the same metric. After division, each cell must grow for a minimum of 5 frames before dividing again.
Please see configuring LAPLinker for more details on these parameters.
-
Insert the code into a for loop.
Now that we have the foundation of the code, we can create a for loop to loop over each frame of the image. Note that the example image has 17 frames.
The full code is posted below. Here, we show a minimal version with comments showing where the different parts of the code would go. However, note the position where
assignToTrack
is called.%Create a LAPLinker object L = LAPLinker; %Set up the tracking parameters L.LinkedBy = 'PixelIdxList'; %...etc. %Process each frame for iT = 1:17 %Read in image I = imread('demo_cell.tif', iT); %Segment the cells mask = mask = I < 1.65e4; %...etc. %Measure data data = regionprops(mask, 'Centroid', 'PixelIdxList', 'MajorAxisLength'); %Assign the data to tracks L = assignToTrack(L, iT, data); end
Visualizing tracking results is one way to validate that the code is operating as expected. The following describes how to generate an output showing both the cell, an outline of the mask, and the object IDs.
-
Compute the perimeter of the mask
maskPerimeter = bwperim(mask);
-
Make the line thicker by dilating the perimeter mask
maskPerimeter = imdilate(maskPerimeter, ones(2));
-
We want to display the outlines in green, so generate an RGB matrix
maskDisplay = cat(3, zeros(size(maskDisplay)), maskDisplay, zeros(size(maskDisplay)));
-
Generate a composite image of both image and mask
Iout = double(imfuse(imadjust(I), maskDisplay, 'blend')); Iout = uint8(Iout ./ max(Iout(:)) * 255);
-
Label the tracks. Since we only want to look at tracks which are in the current frame, we can use the
activeTrackIDs
property.%Insert the IDs of active tracks for iActive = L.activeTrackIDs trackData = getTrack(L, iActive); Iout = insertText(Iout, trackData.Centroid(end, :), trackData.ID, ... 'BoxOpacity', 0, 'TextColor', 'white', 'AnchorPoint', 'Center'); end
-
Show the image
imshow(Iout)
We have now covered all the basic concepts. The following is the full code listing of the example. You can copy and paste the code as a new script.
%Create a LAPLinker object
L = LAPLinker;
%Set up the tracking parameters to use overlapping pixels
L.LinkedBy = 'PixelIdxList';
L.LinkCostMetric = 'pxintersect';
L.LinkScoreRange = [1 12];
L.TrackDivision = true;
L.DivisionParameter = 'PixelIdxList';
L.DivisionScoreMetric = 'pxintersect';
L.DivisionScoreRange = [1 12];
L.MinFramesBetweenDiv = 5;
for iT = 1:17
%Read in image
I = imread('demo_cell.tif', iT);
%Segment the cells
thLvl = 1.65e4;
%Compute mask and fill holes
mask = I < thLvl;
mask = imopen(mask, strel('disk', 3));
mask = imfill(mask, 'holes');
%Separate the cell clumps using watershedding
dd = -bwdist(~mask);
dd(~mask) = -Inf;
dd = imhmin(dd, 1);
LL = watershed(dd);
mask(LL == 0) = 0;
mask = bwmorph(mask, 'thicken', 3);
mask = bwareaopen(mask, 300);
%Measure data
data = regionprops(mask, 'Centroid', 'PixelIdxList', 'MajorAxisLength');
%Assign the data to tracks
L = assignToTrack(L, iT, data);
%Generate an output image to validate
maskDisplay = imdilate(bwperim(mask), ones(2)); %Get outline of the cells
maskDisplay = cat(3, zeros(size(maskDisplay)), maskDisplay, zeros(size(maskDisplay)));
Iout = double(imfuse(imadjust(I), maskDisplay, 'blend'));
Iout = uint8(Iout ./ max(Iout(:)) * 255);
%Insert the IDs of active tracks
for iActive = L.activeTrackIDs
trackData = getTrack(L, iActive);
Iout = insertText(Iout, trackData.Centroid(end, :), trackData.ID, ...
'BoxOpacity', 0, 'TextColor', 'white', 'AnchorPoint', 'Center');
end
%Write the output as a GIF
[imind, cm] = rgb2ind(Iout,256);
if iT == 1
imwrite(imind, cm, 'output.gif', 'Loopcount',inf, 'DelayTime', 0.2);
else
imwrite(imind, cm, 'output.gif', 'writeMode', 'append', 'DelayTime', 0.2);
end
end
This code should produce an output GIF that looks like:
Now that tracking is complete, it is a good idea to save the data.
trackdata = L.tracks;
save('demodata.mat', 'trackdata')