Determining the point-spread-function from a bead image by averaging#

In order to deconvolve a microsocpy image properly, we should determine the point-spread-function (PSF) of the microscope.

See also

import numpy as np
from import imread, imsave
from pyclesperanto_prototype import imshow
import pyclesperanto_prototype as cle
import pandas as pd
import matplotlib.pyplot as plt

The example image data used here was acquired by Bert Nitzsche and Robert Haase (both MPI-CBG at that time) at the Light Microscopy Facility of MPI-CBG. Just for completeness, the voxel size is 0.022x0.022x0.125 µm^3.

bead_image = imread('../../data/Bead_Image1_crop.tif')
(41, 150, 150)

Our example image shows fluorescent beads, ideally with a diameter smaller than the resolution of the imaging setup. Furthermore, the beads should emit light in the same wavelength as the sample we would like to deconvolve later on. In the following image crop we see four fluorescent beads. It is recommended to image a larger field of view, with at least 25 beads. Also make sure that the beads do not stick to each other and are sparsely distributed.

imshow(cle.maximum_x_projection(bead_image), colorbar=True)
imshow(cle.maximum_y_projection(bead_image), colorbar=True)
imshow(cle.maximum_z_projection(bead_image), colorbar=True)
../_images/extract_psf_5_0.png ../_images/extract_psf_5_1.png ../_images/extract_psf_5_2.png

For determining an average PSF, technically we can crop out all the individual beads, align them and then average the images. Therefore, we segment the objects and determine their center of mass.

# Segment objects
label_image = cle.voronoi_otsu_labeling(bead_image)
imshow(label_image, labels=True)
# determine center of mass for each object
stats = cle.statistics_of_labelled_pixels(bead_image, label_image)

df = pd.DataFrame(stats)
df[["mass_center_x", "mass_center_y", "mass_center_z"]]
mass_center_x mass_center_y mass_center_z
0 30.107895 73.028938 23.327475
1 44.293156 111.633430 23.329062
2 76.092850 82.453033 23.299677
3 125.439606 35.972496 23.390951

PSF averaging#

Next, we will iterate over the beads and crop them out by translating them into a smaller PSF image.

# configure size of future PSF image
psf_radius = 20
size = psf_radius * 2 + 1

# initialize PSF
single_psf_image = cle.create([size, size, size])
avg_psf_image = cle.create([size, size, size])

num_psfs = len(df)
for index, row in df.iterrows():
    x = row["mass_center_x"]
    y = row["mass_center_y"]
    z = row["mass_center_z"]
    print("Bead", index, "at position", x, y, z)
    # move PSF in right position in a smaller image
    cle.translate(bead_image, single_psf_image, 
                  translate_x= -x + psf_radius,
                  translate_y= -y + psf_radius,
                  translate_z= -z + psf_radius)

    # visualize
    fig, axs = plt.subplots(1,3)    
    imshow(cle.maximum_x_projection(single_psf_image), plot=axs[0])
    imshow(cle.maximum_y_projection(single_psf_image), plot=axs[1])
    imshow(cle.maximum_z_projection(single_psf_image), plot=axs[2])
    # average
    avg_psf_image = avg_psf_image + single_psf_image / num_psfs
Bead 0 at position 30.107894897460938 73.02893829345703 23.32747459411621
Bead 1 at position 44.293155670166016 111.63343048095703 23.32906150817871
Bead 2 at position 76.09284973144531 82.45303344726562 23.2996768951416
Bead 3 at position 125.43960571289062 35.972496032714844 23.39095115661621
../_images/extract_psf_10_1.png ../_images/extract_psf_10_2.png ../_images/extract_psf_10_3.png ../_images/extract_psf_10_4.png

The average PSF then looks like this:

fig, axs = plt.subplots(1,3)    
imshow(cle.maximum_x_projection(avg_psf_image), plot=axs[0])
imshow(cle.maximum_y_projection(avg_psf_image), plot=axs[1])
imshow(cle.maximum_z_projection(avg_psf_image), plot=axs[2])
avg_psf_image.min(), avg_psf_image.max()
(0.0, 94.5)

After we determined a well centered PSF, we can save it for later re-use. Before doing that, we normalize the PSF. Goal is to have an image where the total intensity is 1. This makes sure that an image that is deconvolvolution using this PSF later on does not modify the image’s intensity range.

normalized_psf = avg_psf_image / np.sum(avg_psf_image)

imshow(normalized_psf, colorbar=True)
normalized_psf.min(), normalized_psf.max()
(0.0, 0.0006259646)
imsave('../../data/psf.tif', normalized_psf)
C:\Users\rober\AppData\Local\Temp\ipykernel_16716\ UserWarning: ../../data/psf.tif is a low contrast image
  imsave('../../data/psf.tif', normalized_psf)