import numpy as np

from IPython.display import Image, display


# Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
%matplotlib inline

import sys
sys.path.append('../') # Add the course directory to the Python path
import cphy.plotting as cplot

Microscopy data and deconvolution methods

Many microscopy techniques are based on the detection of light emitted by fluorescent molecules. The light emitted by a single molecule is very weak and cannot be detected directly. Instead, the sample is illuminated with a strong light source, and the light emitted by the molecules is detected by a camera.

The image obtained in this way is a convolution of the true image with the point spread function (PSF) of the microscope. The PSF is the image of a single molecule under the microscope, and it can often be calculated experimentally or inferred from the optical design imaging setup. Mathematically, the observed image has the form of a convolution between the true image and the PSF:

$$ \begin{align} I_{obs}(x,y) &= (I_{true} * PSF)(x,y) \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} I_{true}(x',y') PSF(x-x',y-y') dx' dy' \end{align} $$

Note the similarity between this equation and the Green's function of a physical system.

In order to increase the effective resolution of the microscope, it is necessary to deconvolve the image with the PSF. This is an ill-posed problem, because there are many different methods for solving it.

fpath = "../resources/fluorescence_image.npy"
im = np.load(fpath, allow_pickle=True)

plt.imshow(im)
cplot.vanish_axes()

## Deconvolve the image

import scipy.fftpack# import fft2, ifft2, fftshift, ifftshift

from scipy.signal import convolve2d

def gaussian_kernel(size: int, sigma: float):
    kernel = np.fromfunction(
        lambda x, y: (1/ (2 * np.pi * sigma**2)) * 
                     np.exp(- ((x - (size - 1) / 2) ** 2 + (y - (size - 1) / 2) ** 2) / (2 * sigma**2)),
        (size, size)
    )
    return kernel / np.sum(kernel)

def richardson_lucy(image, sigma=1, iterations=50):

    image = image.astype(np.float64) #/ np.max(image)

    psf = gaussian_kernel(3, sigma)

    estimate = np.copy(image)
    for i in range(iterations):
        q = image / (1e-8 +convolve2d(estimate, psf, 'same'))
        estimate *= convolve2d(q, psf[::-1, ::-1], 'same')
    return estimate.T


all_channels = []
for i in range(3):
    channel = richardson_lucy(im[..., i])
    channel = 255 * channel / np.max(channel)
    all_channels.append(channel)
restored_image = np.array(all_channels).T
restored_image = restored_image.astype(np.uint8)

# Show image a and restored side-by-side

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(im)
plt.title("Original")
cplot.vanish_axes()

plt.subplot(122)
plt.imshow(restored_image)
plt.title("Deconvolved")
cplot.vanish_axes()

# Compare just green channel to highlight what changed
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(im[..., 1])
plt.title("Original")
plt.subplot(122)
plt.imshow(restored_image[..., 1])
plt.title("Deconvolved")

Text(0.5, 1.0, 'Deconvolved')