Multiscale localisation in 3D and confocal images deconvolution: Howto

Following last week's impulse, here is a tutorial about the subtleties of multiscale particle localisation from confocal 3D images.

For reference, this has been published as

A novel particle tracking method with individual particle size measurement and its application to ordering in glassy hard sphere colloids.

Leocmach, M., & Tanaka, H.

Soft Matter (2013), 9, 1447–1457.

doi:10.1039/c2sm27107a

ArXiV

The corresponding code is available at Sourceforge, including the IPython notebook corresponding to the present tutorial.

Particle localisation in a confocal image¶

First, import necessary library to load and display images

import numpy as np
from scipy.ndimage import gaussian_filter
from matplotlib import pyplot as plt
from colloids import track
%matplotlib inline

We also need to define functions to quickly display tracking results.

def draw_circles(xs, ys, rs, **kwargs):
    for x,y,r in zip(xs,ys,rs):
        circle = plt.Circle((x,y), radius=r, **kwargs)
        plt.gca().add_patch(circle)
        
def display_cuts(imf, centers, X=30, Y=25, Z=30):
    """Draw three orthogonal cuts with corresponding centers"""
    plt.subplot(1,3,1);
    draw_circles(centers[:,0], centers[:,1], centers[:,-2], facecolor='none', edgecolor='g')
    plt.imshow(imf[Z], 'hot',vmin=0,vmax=255);
    plt.subplot(1,3,2);
    draw_circles(centers[:,0], centers[:,2], centers[:,-2], facecolor='none', edgecolor='g')
    plt.imshow(imf[:,Y], 'hot',vmin=0,vmax=255);
    plt.subplot(1,3,3);
    draw_circles(centers[:,1], centers[:,2], centers[:,-2], facecolor='none', edgecolor='g')
    plt.imshow(imf[:,:,X], 'hot',vmin=0,vmax=255);

Raw confocal image¶

Load the test data, which is only a detail of a full 3D stack. Display different cuts.

Note that raw confocal images are heavily distorted in the $z$ direction. We actually have 4 particles in the image: two small on top of each other, one on the side and a slightly larger below and on the side.

imf = np.fromfile('../../../multiscale/test_input/Z_elong.raw', dtype=np.uint8).reshape((55,50,50))
plt.subplot(1,3,1).imshow(imf[30], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,2).imshow(imf[:,25], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,3).imshow(imf[:,:,30], 'hot',vmin=0,vmax=255);

Create a finder object of the same shape as the image. Seamlessly work with 3D images.

finder = track.MultiscaleBlobFinder(imf.shape, Octave0=False, nbOctaves=4)

Feed the image in the tracker with default tracking parameters. The output has four columns : $(x,y,z,r,i)$ where $r$ is a first approximation of the radius of the particle (see below) and $i$ is a measure of the brightness of the particle.

centers = finder(imf)
print centers.shape
print "smallest particle detected has a radius of %g px"%(centers[:,-2].min()) 
display_cuts(imf, centers)

(2, 5)
smallest particle detected has a radius of 4.40672 px

This is a pretty bad job. Only two particle detected correctly. The top and side small particles are not detected.

Let us try not to remove overlapping.

centers = finder(imf, removeOverlap=False)
print centers.shape
print "smallest particle detected has a radius of %g px"%(centers[:,-2].min()) 
display_cuts(imf, centers)

(4, 5)
smallest particle detected has a radius of 4.40672 px

Now we have all the particles.

We have two effects at play here that fool the overlap removal: - The radii guessed by the algorithm are too large. - The particles aligned along $z$ are localised too close together.

Both effects are due to the elongation in $z$.

Simple deconvolution of a confocal image¶

The image $y$ acquired by a microscope can be expressed as

$y = x \star h + \epsilon,$

where $\star$ is the convolution operator, $x$ is the perfect image, $h$ is the point spread function (PSF) of the microscope and $\epsilon$ is the noise independent of both $x$ and $h$. The process of estimating $x$ from $y$ and some theoretical or measured expression of $h$ is called deconvolution. Deconvolution in the presence of noise is a difficult problem. Hopefully here we do not need to reconstruct the original image, but only our first Gaussian blurred version of it. Indeed, after a reasonable amount of blur in three dimensions, the noise can be neglected and we thus obtain:

$y_0 \approx G_0 \star h,$

or in Fourier space

$\mathcal{F}[y_0] = \mathcal{F}[G_0] \times \mathcal{F}[h].$

Once $\mathcal{F}[h]$ is known the deconvolution reduces to a simple division in Fourier space.

Measuring the kernel¶

$\mathcal{F}[h]$ should be measured in an isotropic system and in the same conditions (magnification, number of pixels, indices of refraction, etc.) as the experiment of interest.

Since a 3D image large enough to be isotropic is too big to share in an example, I just use the image of interest. Don't do this at home.

imiso = np.fromfile('../../../multiscale/test_input/Z_elong.raw', dtype=np.uint8).reshape((55,50,50))
kernel = track.get_deconv_kernel(imiso)
plt.plot(kernel);

The kernel variable here is actually $1/\mathcal{F}[h]$

Kernel test¶

What not to do: deconvolve the unblurred image.

imd = track.deconvolve(imf, kernel)
plt.subplot(1,3,1).imshow(imd[30], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,2).imshow(imd[:,25], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,3).imshow(imd[:,:,30], 'hot',vmin=0,vmax=255);

Noisy isn't it? Yet we have saturated the color scale between 0 and 255, otherwise artifacts would be even more apparent.

What the program will actually do: deconvolved the blurred image.

imb = gaussian_filter(imf.astype(float), 1.6)
plt.subplot(1,3,1).imshow(imb[30], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,2).imshow(imb[:,25], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,3).imshow(imb[:,:,30], 'hot',vmin=0,vmax=255);

imbd = track.deconvolve(imb, kernel)
plt.subplot(1,3,1).imshow(imbd[30], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,2).imshow(imbd[:,25], 'hot',vmin=0,vmax=255);
plt.subplot(1,3,3).imshow(imbd[:,:,30], 'hot',vmin=0,vmax=255);

The way the noise is amplified by deconvolution is obvious when we compare the histograms of the 4 versions of the image.

for i, (l, im) in enumerate(zip(['raw', 'blurred', 'deconv', 'blur deconv'], [imf, imb, imd, imbd])):
    plt.plot(
             np.arange(-129,363,8), 
             np.histogram(im, np.arange(-129,371,8))[0]*10**i, 
             label=l+' (%d,%d)'%(im.min(), im.max()));
plt.xlim(-130,360);
plt.yscale('log');
plt.legend(loc='upper right');

Obviously negative values should be discarded, and will be discarded by the algorithm.

Localisation using deconvolved image¶

Not removing the overlaping particles, we use the deconvolution kernel.

reload(track)
centers = finder(imf, removeOverlap=False, deconvKernel=kernel)
print centers.shape
print "smallest particle detected has a radius of %g px"%(centers[:,-2].min()) 
display_cuts(imf, centers)

(5, 5)
smallest particle detected has a radius of 3.79771 px

All particles are found and actually they are not overlaping. However we got a spurious detection, probably some noise that got amplified by the deconvolution. It is easily removed by a threshold in intensity.

print centers

[[ 16.31334384  44.07373885  21.50444345   3.79770755  -0.37019082]
 [ 26.79860906  22.0435863   21.89833819   3.89324706 -15.95536745]
 [ 26.92095001  21.83422163  32.72685299   3.85903411 -17.3271687 ]
 [ 26.11610631  31.30829898  36.11323148   4.33867807 -21.06731591]
 [ 31.77923381  26.78655827  44.05564076   4.03334069 -20.1763175 ]]

centers = centers[centers[:,-1]<-1]
display_cuts(imf, centers)

We can also display the results with respect to the blurred and deconvolved version of the image that we can find by introspection.

display_cuts(finder.octaves[1].layersG[0], centers)

From scales to sizes¶

This is pretty much like in 2D.

In [17]:

from colloids.particles import get_bonds
s = track.radius2sigma(centers[:,-2], dim=3)
bonds, dists = get_bonds(positions=centers[:,:-2], radii=centers[:,-2], maxdist=3.0)
brights1 = track.solve_intensities(s, bonds, dists, centers[:,-1])
radii1 = track.global_rescale_intensity(s, bonds, dists, brights1)
centers2 = np.copy(centers)
centers2[:,-2] = radii1
display_cuts(imf, centers2)