PGURE-SVT Demonstration

Tom Furnival (tjof2@cam.ac.uk)

PGURE-SVT is an algorithm designed to denoise image sequences acquired in microscopy. It exploits the correlations between consecutive frames to form low-rank matrices, which are then recovered using a technique known as nuclear norm minimization. An unbiased risk estimator for mixed Poisson-Gaussian noise is used to automate the selection of the regularization parameter, while robust noise and motion estimation maintain broad applicability to many different types of microscopy.

You can read more about the algorithm and its applications in:

T Furnival, RK Leary, PA Midgley "Denoising time-resolved microscopy sequences with singular value thresholding." (Article in press). DOI:10.1016/j.ultramic.2016.05.005

The source code and compiled Windows libraries are available to download from http://tjof2.github.io/pgure-svt/.

Use with HyperSpy

This example notebook shows how PGURE-SVT can be combined with HyperSpy, which is an open-source Python library that makes signal handling and processing straightforward in Python, with a friendly API. While you can use pguresvt.pguresvt.SVT to denoise a NumPy array directly, pguresvt.hspysvt.HSPYSVT can instead denoise a HyperSpy signal.



In [ ]:

    
# Configures the plotting backend
%matplotlib qt

# Import NumPy and the HyperSpy API
import numpy as np
import hyperspy.api as hs

# Import the HyperSpy wrapper for PGURE-SVT
from pguresvt import hspysvt
from pguresvt.pguresvt import PoissonGaussianNoiseGenerator

1. Simulated dataset

First, we load the simulated dataset using HyperSpy.



In [ ]:

    
# Load example dataset
movie = hs.load("examplesequence.tif")

# Truncate to 25 frames, and plot the result
movie = movie.inav[:25]
movie.plot(navigator='slider')

Now we corrupt the movie with using a noise generator for mixed Poisson-Gaussian noise, according to the equation: where the parameters are defined as:

alpha = detector gain
mu = detector offset
sigma = detector noise



In [ ]:

    
# Detector gain
detector_gain = 0.1
# Detector offset
detector_offset = 0.1
# Detector noise
detector_sigma = 0.1

noisy = PoissonGaussianNoiseGenerator(movie._data_aligned_with_axes,
                                      alpha=detector_gain,
                                      mu=detector_offset,
                                      sigma=detector_sigma)
noisy_movie = hs.signals.Image(noisy)
noisy_movie.plot(navigator='slider')

Next we initialise the SVT denoising function. The full list of options (with default values) is:

hspysvt.HSPYSVT(patchsize=4,
                patchoverlap=1,
                length=15,
                optimize=True,
                threshold=0.5,
                estimatenoise=True,
                alpha=-1., 
                mu=-1., 
                sigma=-1.,                
                arpssize=7, 
                tol=1e-7,
                median=5,
                hotpixelthreshold=10,
                numthreads=1)

In this example we do not use the noise estimation procedure, and instead provide the known parameters to the algorithm directly. This information is used by the PGURE optimizer to calculate the threshold.

Note: If you have a multicore machine, you can set numthreads > 1 to speed up the calculation.



In [ ]:

    
# Initialize with suggested parameters
svt = hspysvt.HSPYSVT(patchsize=4,
                      estimatenoise=False,
                      alpha=detector_gain,
                      mu=detector_offset,
                      sigma=detector_sigma,
                      tol=1e-5,
                      numthreads=2)

Now we are able to run the denoising and plot the result:



In [ ]:

    
# Run the denoising
denoised_movie = svt.denoise(noisy_movie)

# Plot denoised data
denoised_movie.plot(navigator='slider')

2. Time-resolved ADF-STEM image sequence

In this example we apply PGURE-SVT to an experimental dataset of a nanoparticle acquired using ADF-STEM. This image sequence contains 51 frames at a rate of 4 frames per second. The results of this denoising are shown in Fig. 11 of the paper.

For larger images, such as the 256x256 pixels here, you can use the patchoverlap parameter to control the trade-off between speed and accuracy of the denoising procedure. This reduces the number of patches the algorithm works with, at the expense of introducing possible edge artefacts between patches.

For the experimental sequence, the detector offset (mu) was known beforehand, so a noise estimation procedure is used for the other values.



In [ ]:

    
# Load example dataset and plot
expt_movie = hs.load("experimentalsequence-NP.tif)
expt_movie.plot(navigator='slider')

# Initialize with suggested parameters, optimized for speed
expt_svt = hspysvt.HSPYSVT(patchsize=4,
                           patchoverlap=2,
                           mu=0.075,
                           numthreads=2)

# Run the denoising
denoised_movie = expt_svt.denoise(expt_movie)

# Plot denoised data
denoised_movie.plot(navigator='slider')

The parameters used to generate Fig. 11 in the paper are below. Note that using these values can be slow, taking ~30 seconds per frame.

expt_svt = hspysvt.HSPYSVT(patchsize=4,
                           patchoverlap=1,
                           mu=0.075,
                           tol=1e-8,
                           arpssize=11,
                           numthreads=4)