Feature finding in 3D confocal images¶

Setup IPython, pims, trackpy¶

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('figure',  figsize=(10, 6))
import numpy as np
import pandas as pd
from pandas import DataFrame, Series  # for convenience

import pims
import trackpy as tp

Import 3D images¶

We make use of ImageSequenceND. This reader class searches for files matching a certain search string. The filename should contain indentifiers for each axis (dimension). For example "demo_t001_z001.png" contains identifiers t and z. We set the returned image axes to zyx and the iterating axis to t. In this way we get 3D images indexed by t.

For this tutorial, the images are inside a zipfile. You can extract the images from the specified archive to check them out.

frames = pims.ImageSequenceND(r'../sample_data/pmma_colloids.zip', axes_identifiers = ['z', 't'])
frames.bundle_axes = ['z', 'y', 'x']
frames.iter_axes = 't'
frames

<ImageSequenceND>
Source: D:GebruikersCasperDocumentenGitHubtrackpy-examplessample_datapmma_colloids.zip
Axes: 4
Axis 'y' size: 128
Axis 'x' size: 128
Axis 'z' size: 25
Axis 't' size: 33
Pixel Datatype: uint8

The time-index of the first frame is always 0, regardless of the index in the filename. We can use pims 3d scrollable stack feature to display a 3D image. Use the scroll wheel to scroll through the focal planes.

This feature will not work online, download the example and run it yourself to checkout the scrollable stack feature.

frames[0]

Feature finding¶

Now we can use trackpy.locate to generate a table (a pandas DataFrame) of the feature coordinates. locate only requires a feature diameter. This needs to be an odd integer and can be specified for each dimension seperately. The order of dimensions is the same as the order of the image array, typically (z, y, x).

features = tp.locate(frames[0], diameter=(5, 9, 9))
features.head()  # displays first 5 rows

Because we specified different diameters for each dimension, trackpy we return size and ep (static error) for each dimension seperately. We can plot the feature locations on top of the original picture. You can see that there are no features found at the edges. This is because these particles partly fall out of the measurement box.

tp.annotate3d(features, frames[0])

To see wether the features are biased by the pixelation, we can take histograms of the decimal part of the coordinates. As particles don’t know about the pixel grid, these distributions should be flat.

tp.subpx_bias(features)
print 'Features found: {0}'.format(len(features))

C:Anacondalibsite-packagespandastoolsplotting.py:2954: UserWarning: To output multiple subplots, the figure containing the passed axes is being cleared
  "is being cleared", UserWarning)

Features found: 144

We can do better by increasing the diameter size:

features2 = tp.locate(frames[0], diameter=(7, 11, 11))
tp.subpx_bias(features2)
print 'Features found: {0}'.format(len(features2))

Features found: 87

But now we lost 54 features! This is because of two things: 1) the margins become larger, 2) trackpy removes particles that are closer than diameter from each other. Mainly the second effect is responsible for the large loss in this case. We can force the minimum separation in trackpy to a lower value, to keep features that are closer than diameter from each other. We should keep in mind that the coordinates of two close particles will be biased towards each other.

features3 = tp.locate(frames[0], diameter=(7, 11, 11), separation=(3, 11, 11))
print 'Features found: {0}'.format(len(features3))
tp.annotate3d(features3, frames[0])

Features found: 165

Batch processing¶

We use the optimized settings and iterate through the full movie.

f = tp.batch(frames, diameter=(7, 11, 11), separation=(3, 11, 11))

Frame 32: 159 features

Linking¶

The linking features of trackpy support an arbitrary number of dimensions. The argument pos_columns will determine which and how many dimensions will be used. The standard set of columns is ['x', 'y']. For 3d we have to set this to ['x', 'y', 'z']. Linking treats all dimensions equally. For this example we are using a file with unequal pixel sizes: 0.21 microns in XY and 0.75 microns in Z, so we first define new columns in units of microns.

f['xum'] = f['x'] * 0.21
f['yum'] = f['y'] * 0.21
f['zum'] = f['z'] * 0.75

Experimentally, there is a tradeoff between framerate and z pixel size. If you want to have higher z resolution, you scan slower and your framerate goes down. In this specific experiment the z resolution is low, so that we are able to link particles between subsequent frames. We try several linking distances here.

for search_range in [1.0, 1.5, 2.0, 2.5]:
    linked = tp.link_df(f, search_range, pos_columns=['xum', 'yum', 'zum'])
    hist, bins = np.histogram(np.bincount(linked.particle.astype(int)),
                              bins=np.arange(30), normed=True)
    plt.step(bins[1:], hist, label='range = {} microns'.format(search_range))
plt.gca().set(ylabel='relative frequency', xlabel='track length (frames)')
plt.legend()
pass

Frame 32: 159 trajectories present

Mean squared displacement¶

We use a search range of 2.0 microns, because there is no significant improvement between 2.0 and 2.5 microns.

linked = tp.link_df(f, 2.0, pos_columns=['xum', 'yum', 'zum'])

Frame 32: 159 trajectories present

With the trajectories we can calculate the mean squared displacement. To calculate the mean squared displacement from all three dimensions, we need to tell trackpy which position columns to use. It defaults to ['x', 'y'].

We specify mpp=1 because we have already transformed our data to microns.

msd3D = tp.emsd(linked, mpp=1, fps=1/0.8582, max_lagtime=20,
                pos_columns=['xum', 'yum', 'zum'])
ax = msd3D.plot(style='o', label='MSD in 3D')
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]', xlabel='lag time $t$')
ax.set(xlim=(0, 16), ylim=(0, 20))
ax.legend(loc='upper left')
pass

The mean squared displacement in three dimensional Brownian motion is described by the following relation:

Let’s do a fit and calculate the diffusion constant!

slope = np.linalg.lstsq(msd3D.index[:, np.newaxis], msd3D)[0][0]

ax = msd3D.plot(style='o', label='MSD in 3D')
ax.plot(np.arange(20), slope * np.arange(20), label='linear fit')
ax.set(ylabel=r'$\langle \Delta r^2 \rangle$ [$\mu$m$^2$]', xlabel='lag time $t$')
ax.set(xlim=(0, 16), ylim=(0, 20))
ax.legend(loc='upper left')
print(r'The diffusion constant is {0:.2f} μm²/s'.format(slope / 6))

The diffusion constant is 0.18 μm²/s

About this work¶

This tutorial was written by Casper van der Wel, as part of his PhD thesis in Daniela Kraft’s group at the Huygens-Kamerlingh-Onnes laboratory, Institute of Physics, Leiden University, The Netherlands. This work was supported by the Netherlands Organisation for Scientific Research (NWO/OCW).