Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Variance lower bound on fluorescence microscopy image denoising

Open Access Open Access

Abstract

The signal to noise ratio of high-speed fluorescence microscopy is heavily influenced by photon counting noise and sensor noise due to the expected low photon budget. Denoising algorithms are developed to decrease these noise fluctuations in microscopy data by incorporating additional knowledge or assumptions about imaging systems or biological specimens. One question arises: whether there exists a theoretical precision limit for the performance of a microscopy denoising algorithm. In this paper, combining Cramér-Rao Lower Bound with constraints and the low-pass-filter property of microscope systems, we develop a method to calculate a theoretical variance lower bound of microscopy image denoising. We show that this lower bound is influenced by photon count, readout noise, detection wavelength, effective pixel size and the numerical aperture of the microscope system. We demonstrate our development by comparing multiple state-of-the-art denoising algorithms to this bound. This method establishes a framework to generate theoretical performance limit, under a specific prior knowledge, or assumption, as a reference benchmark for microscopy denoising algorithms.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Rapid development of fluorescence microscopy techniques together with the availability of fast and sensitive cameras are enabling molecular observations at unprecedented spatial and temporal resolution. High-speed imaging of fluorescently tagged molecules suffers from low signal to noise ratio (SNR) due to the limited photon budget per fluorescent emitter. When ignoring sensor noise, the detected photons in each pixel follow Poisson distribution [1]. As a result, the SNR of a microscopy image decreases rapidly with the decreasing number of detected photons. At these low light conditions, quantitative biological measurements result in imprecisions — fluctuations of measured signal mainly come from the uncertainty of photon detection rather than the underlying biological processes.

To improve the measurement precision, denoising algorithms are developed to decrease the noise fluctuation in the microscopy data caused by the photon detection process and the sensor. The key to such noise reducing capacity is the incorporation of additional prior knowledge or assumptions about imaging systems or biological specimens [27], for example, through exploiting the low pass filter property of a well-designed microscope [3] or assuming self-similarities of local structures in the specimen [47].

An ideal denoising algorithms will maintain the quantitative nature of microscopy images by providing estimations of pixel intensities with low variance and minimum biases or artifacts. However, common among the denoising algorithms which we tested, there are tradeoffs between the estimation variance and the introduced bias. Typically, one could shift denoising performance by increasing bias in trade for decreasing the denoising variance. Regarding this tradeoff, one basic question arises: How precise is precise enough? Is there a precision limit which a denoising algorithm could achieve at its best given a particular prior knowledge and a certain bias level?

Previously, Levin et al. [8] proposed a Bayesian approach to calculate the mean square error lower bound of natural image denoising algorithms. Without formulating a parametric constraint, the prior knowledge (or distribution) of image structures is drawn from a collection of 1010 natural images. This approach results in a useful tool to calculate overall performance limit of a denoising algorithm based on a generic natural image prior. However, microscopy images vary significantly with specimen types, magnifications, pixel sizes, objective types as well as sensor types. There is no universal prior that can take both the specimen structure and the imaging system into account. This calls for a prior that is generally valid for images acquired in microscopes, independent of specimen structures.

In this context, we developed a mathematical framework to calculate the theoretical variance lower bound of unbiased denoising algorithms for microscopy images. Specifically, we provided an example of this lower bound based on the low-pass-filter property of microscopes serving as a reference benchmark for multiple microscopy denoising algorithms. This bound also allows us to explore the expected performance limit of denoising tasks in relationship with the detected photons (e.g. wavelength ($\lambda $) and expected photon count) and the microscope system (e.g. objective’s numerical aperture ($NA$), and effective camera pixel size) given an unbiased denoising algorithm.

2. Methods

2.1. Notation

To simplify notation, throughout this paper we used bold upper case letter to denote matrix (e.g. ${\boldsymbol{A}}$), bold lower case letter to denote vector (e.g. ${\boldsymbol a}$) and normal letter despite of upper or lower case to denote scalar (e.g. a or $A$). All variables were denoted in italic fonts while operations were denoted in normal fonts (e.g. ${\mathbf A}$).

2.2. Frequency constraint of microscopy images

Our aim is to develop a variance lower bound for unbiased microscopy denoising algorithms by considering a common property of a far-field optical microscope system: its frequency response is characterized by the optical transfer function (OTF) [9]. In a typical microscope system, the OTF boundary (with the radius of $2NA/\lambda $) defines the region of detectable spatial frequencies whereas frequencies outside this boundary cannot transmit through the microscope. Mathematically, this low-pass-filter property of microscopy systems forms a rather strict constraint: The 2D Fourier transform of the underlying transmitted signal through a microscope must be zero outside the OTF boundary (Fig. 1).

 figure: Fig. 1.

Fig. 1. Schematic diagram of obtaining cCRLB, bias, and variance maps for evaluating performance of microscopy image denoising algorithms. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. CRLB is the variance lower bound on estimation of ideal image calculated from the noise model [Eq. (4)–(10)]. OTF is the optical transfer function. cCRLB is the variance lower bound on estimation of ideal image while taking CRLB and additional prior knowledge on frequency constraint [Eq. (3)] into consideration. Bias map and variance map (on the right) were calculated pixel by pixel from 100 frames of denoised images (using NCS denoising algorithm).

Download Full Size | PDF

To express this frequency constraint analytically, we use discrete Fourier transform (DFT) matrix to express the 2D Fourier transform of the image. For each frequency component of the image, we have

$$\begin{aligned}{[{FT\{{\boldsymbol U} \}} ]_{p,q}} & = \mathop \sum \nolimits_{j = 0}^{n - 1} \mathop \sum \nolimits_{k = 0}^{n - 1} {U_{j,k}}{e^{ - i\frac{{2\pi }}{n}({pj + qk} )}} \\ & = \mathop \sum \nolimits_{j = 0}^{n - 1} \mathop \sum \nolimits_{k = 0}^{n - 1} {W_{p,j}}{U_{j,k}}{W_{k,q}}\\ &= {[{{\boldsymbol{WUW}}} ]_{p,q}} \\ & = 1_{\boldsymbol p}^\textrm{T}{\boldsymbol{WUW}}{1_{\boldsymbol q}} \end{aligned}$$
where $({p,q} )$ is a coordinate in Fourier space of the image. ${\boldsymbol U}$ is the $n\ast n$ 2D array of the ideal image where each element ${U_{j,k}}$ is the expected photon count of the pixel at position $({j,k} )$. $\; FT\{{\boldsymbol U} \}$ is the 2D Fourier transform of ${\boldsymbol U}$. ${1_{\boldsymbol p}}$ is a column vector of size $n$ with ${p^{th}}$ element being equal to 1 while the other elements being equal to 0. Superscript $\textrm{T}$ denotes the conjugate transpose of a matrix. ${\boldsymbol W}$ is the $n\ast n$ DFT matrix defined in [10],
$${W_{p,q}} = {e^{ - i\frac{{2\pi pq}}{n}}}. $$

Thus we express the constraint by setting the amplitude square of components outside the OTF boundary to be 0,

$$\begin{aligned} {E_{p,q}} &= {({1_{\boldsymbol p}^\textrm{T}{\boldsymbol{WUW}}{1_{\boldsymbol q}}} )^\textrm{T}}({1_{\boldsymbol p}^\textrm{T}{\boldsymbol{WUW}}{1_{\boldsymbol q}}} ) \\ & = 1_{\boldsymbol q}^\textrm{T}{{\boldsymbol W}^\textrm{T}}{{\boldsymbol U}^\textrm{T}}{{\boldsymbol W}^\textrm{T}}{1_{\boldsymbol p}}1_{\boldsymbol p}^\textrm{T}{\boldsymbol{WUW}}{1_{\boldsymbol q}} \\ & = 0 \end{aligned}$$
for all pixel coordinates $({p,q} )$ outside the OTF boundary (Fig. 1). Here scalar ${E_{p,q}}$ is the magnitude square of the 2D Fourier transform of the image at $({p,q} )$.

2.3. Noise model of the image

Image denoising can be regarded as an estimation process during which one seek to estimate the noise-free ideal image from a detected noisy image with the help of additional knowledge. We consider the expected photon counts at different pixels as parameters to be estimated,

$${\boldsymbol \mu } = {[{{\mu_1},{\; }{\mu_2},{\; } \ldots ,{\; }{\mu_N}} ]^\textrm{T}}$$
where ${\boldsymbol \mu } = \textrm{vec}({\boldsymbol U} )$ is the vectorization form of ideal image ${\boldsymbol U}$. $N = {n^2}$ is the total number of pixels within a 2D image. Here we also define the observation of the image,
$${\boldsymbol z} = {[{{z_1},{\; }{z_2},{\; } \ldots ,{\; }{z_N}} ]^\textrm{T}}$$
where ${z_j}$ is the observed readout of the ${j^{th}}$ pixel in the noisy image.

To consider photon-counting noise and sensor noise, we assume noise from each pixel is independent and can be modeled as the combination of two noise types [11]:

$${Z_j} = {g_j}{S_j} + {Y_j} + {o_j}$$
where ${Z_j}$, ${S_j}$ and ${Y_j}$ are random variables representing ${j^{th}}$ pixel’s camera readout in unit of analog to digit unit (ADU), detected photoelectrons in unit of e- (hereafter referred as photons) and the pixel-dependent readout noise in unit of ADU. Gain factor ${g_j}$ represents the pixel-dependent gain (unit: ADU/e-) and ${o_j}$ is the pixel-dependent offset pre-engineered into the readout process in order to prevent negative counts from camera sensors. In this work, we assume ${S_j}$ as a Poisson random variable with an expectation of ${\mu _j}$ and ${Y_j}$ as a zero-mean Gaussian random variable with a variance of $\sigma _j^2$.

As a result, the probability of obtaining a specific camera readout ${z_j}$ given an expected photon counts ${\mu _j}$ can be expressed as [11,12],

$$\begin{aligned} {p_{{Z_j}}}({{z_j}|{\mu_j}} ) & = \mathop \sum \nolimits_{\tau = 0}^\infty {p_{{S_j}}}({\tau |{\mu_j}} ){p_{{Y_j}}}({{z_j} - {g_j}\tau - {o_j}} ) \\ & = \mathop \sum \nolimits_{\tau = 0}^\infty \frac{{{e^{ - {\mu _j}}}\mu _j^\tau }}{{\tau !}}\frac{1}{{\sqrt {2\pi } {\sigma _j}}}\textrm{exp}\left\{ { - \frac{{{{({{z_j} - {g_j}\tau - {o_j}} )}^2}}}{{2\sigma_j^2}}} \right\} \end{aligned}$$
where ${p_{{Z_j}}}()$ is the probability density function (PDF) of ${Z_j}$, ${p_{{S_j}}}()$ is the probability mass function of ${S_j}$ and ${p_{{Y_j}}}()$ is the PDF of ${Y_j}$.

2.4. Variance lower bound on denoising using OTF prior

Cramér-Rao Lower Bound [13] (CRLB) characterizes the lower bound on variance of unbiased estimators for parameters based on the likelihood function. In presence of constraints imposed on the parameters to be estimated, the lower bound on their variances can be further reduced since these constraints effectively reduce the dimensions of the parameter space [14]. This constraint-considering lower bound is referred as the constrained Cramér-Rao Lower Bound (cCRLB).

We define the covariance matrix of the parameter vector,

$${\boldsymbol \varSigma }({\boldsymbol \mu } )= {\mathbf E}[{{{({\hat{{\boldsymbol \mu }} - {\boldsymbol \mu }} )}^\textrm{T}}({\hat{{\boldsymbol \mu }} - {\boldsymbol \mu }} )} ]$$
where ${\mathbf E}$ is taking the expectation over the distribution of observations ${\boldsymbol z}$.

According to Cramér and Rao et al. [15,16], the covariance matrix of the estimator satisfies the following inequation when regularity condition is satisfied,

$${\boldsymbol \varSigma }({\boldsymbol \mu } )\ge {{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$$
where matrix inequality means that ${\boldsymbol \varSigma }({\boldsymbol \mu } )- {{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$ is positive semidefinite. ${\boldsymbol I}({\boldsymbol \mu } )$ is the $N \times N$ Fisher information matrix.

In our case, as we assume the noise statistics at different pixels are independent, the off-diagonal elements of ${\boldsymbol I}({\boldsymbol \mu } )$ are zeros while diagonal elements of $\; {\boldsymbol I}({\boldsymbol \mu } )$ can be calculated as [12],

$$[\boldsymbol{I}(\boldsymbol{\mu})]_{j, j}=\frac{1}{\mu_{j}^{2}} \int_{-\infty}^{\infty} \frac{\left[\sum_{\tau=0}^{\infty} \tau p_{S_{j}}\left(\tau \mid \mu_{j}\right) p_{Y_{j}}\left(z_{j}-g_{j} \tau-o_{j}\right)\right]^{2}}{\sum_{\tau=0}^{\infty} p_{S_{j}}\left(\tau \mid \mu_{j}\right) p_{Y_{j}}\left(z_{j}-g_{j} \tau-o_{j}\right)} d z_{j}-1$$
When we take the constraints into account, the cCRLB can be calculated as [14],
$${\boldsymbol \varSigma }({\boldsymbol \mu } )\ge {\boldsymbol P}({\boldsymbol \mu } ){{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$$
${\boldsymbol P}({\boldsymbol \mu } )$ is the projection matrix which takes the following form [14],
$${\boldsymbol P}({\boldsymbol \mu } )= {\boldsymbol I} - {{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )\nabla {{\boldsymbol G}^\textrm{T}}({\boldsymbol \mu } ){[{\nabla {\boldsymbol G}({\boldsymbol \mu } ){{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )\nabla {{\boldsymbol G}^\textrm{T}}({\boldsymbol \mu } )} ]^{ - 1}}\nabla {\boldsymbol G}({\boldsymbol \mu } ). $$
${\boldsymbol I}$ is an identity matrix of size $N\ast N$. ${\boldsymbol G}({\boldsymbol \mu } )$ is a column set of constraint functions which satisfies:
$${\boldsymbol G}({\boldsymbol \mu } )= {[{{G_1}({\boldsymbol \mu } ),{\; }{G_2}({\boldsymbol \mu } ),{\; } \ldots ,{\; }{G_K}({\boldsymbol \mu } )} ]^\textrm{T}} = 0$$
where ${G_j}({\boldsymbol \mu } )$ is a scalar function of ${\boldsymbol \mu }$ and K is the total number of constraints. $0$ refers to a column vector of zeros of size K.

We note the convention that

$${[{\nabla {\boldsymbol G}({\boldsymbol \mu } )} ]_{j,k}} = \frac{{\partial {G_j}({\boldsymbol \mu } )}}{{\partial {\mu _k}}}.$$
Using the OTF constraint shown in Eq. (3), the constraint function can be expressed as ${G_j}({\boldsymbol \mu } )= {E_{p,q}} = 0$ where pixel $({p,q} )$ in the Fourier domain is the ${j^{th}}$ pixel when we arranged all pixels outside OTF boundary in a raster scan order. Accordingly, since ${\boldsymbol \mu } = \textrm{vec}({\boldsymbol U} )$, each row of $\nabla {\boldsymbol G}({\boldsymbol \mu } )$ can be calculated as,
$$\nabla {G_j}({\boldsymbol \mu } )= \textrm{vec}{\left( {\frac{{\partial {E_{p,q}}}}{{\partial {\boldsymbol U}}}} \right)^\textrm{T}}.$$
Here, we note the convention,
$${\left[ {\frac{{\partial {E_{p,q}}}}{{\partial {\boldsymbol U}}}} \right]_{j,k}} = \frac{{\partial {E_{p,q}}}}{{\partial {U_{j,k}}}}.$$
To simplify expression, we define
$${\boldsymbol D} = {{\boldsymbol W}^\textrm{T}}{1_{\boldsymbol p}}1_{\boldsymbol p}^\textrm{T}{\boldsymbol W},\;{\boldsymbol b} = {\boldsymbol W}{1_{\boldsymbol q}},\;{\boldsymbol D^{\prime}} = {\boldsymbol b}{{\boldsymbol b}^{\mathbf T}} = {\boldsymbol W}{1_{\boldsymbol q}}1_{\boldsymbol q}^\textrm{T}{{\boldsymbol W}^\textrm{T}}.$$

Then Eq. (3) can be expressed as,

$$E_{p, q}=\boldsymbol{b}^{\mathrm{T}} \boldsymbol{U}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{U} \boldsymbol{b}$$
As a result, we have [17],
$$\frac{{\partial {E_{p,q}}}}{{\partial {\boldsymbol U}}} = 2{\boldsymbol{DUb}}{{\boldsymbol b}^\textrm{T}} = 2{\boldsymbol{DUD}^{\prime}}$$
In summary, we first calculate the unconstrained CRLB ${{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$ and then the projection matrix ${\boldsymbol P}({\boldsymbol \mu } )$. Finally, the cCRLB is obtained by multiplying the resulting projection matrix with the unconstrained CRLB. We refer our readers to the supplementary software, Code 1, Ref. [18].

We would like to note that the OTF prior represents the common and minimum knowledge available for microscopy image denoising. Algorithms that take advantage of other prior knowledge or assumptions may achieve lower variance than the OTF-based cCRLB.

3. Results and analysis

3.1. Comparison of denoising algorithms’ performance with cCRLB

Here, we compared cCRLB with the achieved performance of three denoising algorithms: noise correction for scientific Complementary Metal–Oxide–Semiconductor (sCMOS) camera (NCS) [3], Noise2Void (N2V) [2] and automatic correction of sCMOS-related noise (ACsN) [5]. These algorithms rely on different principles: noise separations based on OTF, image denoising using neural network and the combination of local structure similarity and noise level estimation. To visualize and quantify the effectiveness of tested denoising algorithms, we evaluated their performances at a low signal-to-noise ratio (SNR) condition with peak expected photon count of 20 per pixel and a uniform background of 10 photons per pixel. The gain map and variance map of CMOS sensor were obtained from an sCMOS camera (Orca-Flash4.0v3, Hamamastu, Japan) shown in (Fig. 9 in the Appendix) through pixel-dependent noise calibration [11]. We quantified denoising performance by calculating the achieved estimation variance as intensity fluctuation of the denoised image sequence and bias as the deviation between the ideal image and the mean of the denoised image sequence. We found that all three tested algorithms reduced the estimated variance substantially compared to calibrated noisy image (Fig. 2). Examining both bias and variance of the denoising algorithms, we found N2V, which was trained with Siemens star structures, resulted in the largest variance improvement (Fig. 2) achieving a mean variance of 1.09 (e-)2 but at the cost of high absolute biases (∼0.68 e- on average at structure region defined as the region above average photon counts). These highly biased estimations cause the denoised image being statistically distorted from the ideal image. ACsN shows significantly reduced bias while achieving a mean variance of 11.34 (e-)2. NCS denoising algorithm (inherent balancing parameter $\alpha = 3$) achieves a similar level of bias as ACsN with a mean variance of 8.87 (e-)2. To demonstrate the precision improvement using denoising algorithms relative to the cCRLB, we calculated the pixel wise ratio between achieved variance of a certain denoising algorithm and cCRLB. We found the resulting variance/cCRLB ratios of NCS were 1.56 on average. The ratios of N2V are 0.12 on average at the background regions while being 0.24 on average at regions with structures. ACsN had the ratio of 1.63 on average at the background region while being 2.24 on average at regions with structures.

 figure: Fig. 2.

Fig. 2. Comparison of three denoising algorithms and cCRLB. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy image (raw noisy image subtracted by the offset and divided by the gain map) and denoised images from NCS ($\alpha = 3$), N2V and ACsN algorithms. Variance map (second row) and bias map (third row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps are shown in the fourth row.

Download Full Size | PDF

It’s worth noting that N2V as a neural network based denoising algorithm does not require ideal images as reference in the training process but solely requires similar noisy images as its training set. When the network was trained based on types of structures that are different from the one we intend to denoise, the performance of denoising varied significantly (Fig. 7 in the Appendix). We tested three types of structures in our training sets: Siemens stars, random block patterns and worm-like-chain model simulated microtubules. In order to cover a wide range of SNR in the training set, we set the intensities of the training sets such that the peak expected photon counts per pixel range from 10 to 40 for all 3 types of structures whereas test images have a peak expected photon count of 20. We found N2V trained with Siemen star structures performed the best with a mean variance of 1.41 (e-)2 and bias of 0.79 (e-) compared to mean variances of 0.96/1.16 (e-)2 and biases of 1.88/1.44 (e-) when trained with random pattern/microtubules. Visually, N2V tended to create similar patterns in the denoised images as the patterns in the training images, and therefore when training and testing structures match, N2V gave the best performance during our tests (Fig. 7 in the Appendix).

To compare the three denoising algorithms at higher photon count conditions (Fig. 8 in the Appendix), we increased the peak expected photon counts to be 50 per pixel, with a uniform background of expected photon counts of 10 per pixel. In this setting, the performance of NCS and N2V (trained with Siemens star structures) were similar with their bias and variance increased accordingly with higher photon counts. However, ACsN’s performance improved significantly as the ratio between its variance and cCRLB achieved to 0.97 on average at structure regions and 0.51 on average at the background regions, with relative low absolute bias (∼ 0.40 e- on average at structure regions). A possible reason of ACsN’s improved denoising performance on variance is that ACsN recognizes the local patterns in the image and used them to denoise similar local patches. However, this assumption on structure similarity wasn’t considered in our cCRLB calculation. At low-photon-counts condition (as shown in Fig. 2), the SNR was low such that ACsN may falsely recognize noisy patches and use them to estimate other patches translating the noise-induced pattern to other image locations which resulted in estimation imprecisions.

3.2. Sensor noise influence on cCRLB

Other than photon counting noise, sensor noise (e.g. pixel-dependent readout noise of a sCMOS sensor) also affects the performance of denoising algorithm. Interestingly, by quantifying the influence of sensor noise on denoising variance lower bound, we found that cCRLB of a single pixel was insensitive towards its intrinsic readout-noise level. This behavior differed from CRLB which increased rapidly with increasing readout-noise standard deviation (SD). To illustrate this observation, we compared cCRLB and CRLB of individual pixel by changing one pixel’s readout noise while keeping the remaining pixels’ readout noise levels unchanged (Fig. 3, pixels at different locations represents independent experiments). As shown in Fig. 3(b), when the pixel had ignorable readout-noise SD (σ=0) and constant gain of 2.17 (Fig. 3(b) Location 1), the cCRLB was as low as 7.66 (e-)2, 3.5 times lower than the CRLB (26.57 (e-)2). When increasing the readout-noise SD of the pixel (σ=40 ADU), the cCRLB increased to 10.53 (e-)2 compared to CRLB 491.3 (e-)2, a 47 times difference.

 figure: Fig. 3.

Fig. 3. Comparison of CRLB and cCRLB at different readout noise levels. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm (a) The ideal image with the target pixels marked. (b) Comparisons between CRLB and cCRLB with respect to increasing readout-noise SD in the target pixel while those of the remaining pixels were kept constant. (c) The ratios of CRLB/cCRLB of the target pixel and the mean ratios of its neighbors (8 surrounding pixels) at different readout noise levels of the target pixel.

Download Full Size | PDF

The ratio between CRLB and cCRLB (Fig. 3(c)) shows that when readout-noise of a pixel increases, the benefit of denoising algorithm for the specific pixel increases as well. Nevertheless, the increased readout noise still lower the information content within the image, and such information loss also deteriorates the achievable precision of denoising at surrounding pixels reflected by an increase of cCRLB of them. We observed these increases of cCRLB were more significant on neighboring pixels than pixels further away. To demonstrate this, we calculated the ratio between CRLB and cCRLB for both the target pixel (which we tuned the readout noise level) and the neighboring pixels (whose readout noise level were kept constant). Figure 3(c) shows that with readout-noise SD of target pixel increased from 0 to 40 (ADU), the mean ratio between CRLB and cCRLB of 8 neighboring pixels decreases from 3.2 to 2.9, a small cost due to the pixel correlation introduced by the constraints (Fig. 10 in the Appendix) paid for the significant improvement of estimation variance at the target pixel (CRLB/cCRLB increases from 2.17 to 47).

This result shows the possibility of denoising algorithm estimating pixels precisely with high readout noise. This further highlights the importance of denoising algorithm development in microscopy images, especially for sCMOS sensor where pixel-dependent readout-noise SD varies significantly (e.g. from 1.5 to 40 (ADU)) among pixels.

3.3. Numerical aperture influence on the denoising lower bound

Further, we demonstrated that cCRLB is heavily influenced by the numerical aperture ($NA$) as well as the effective camera pixel size in the specimen. We expressed the radius of the OTF boundary in the unit of numbers of pixels in the Fourier space as

$${R_{\textrm{OTF}}} = \frac{{2NA}}{\lambda } \cdot {n_{\textrm{img}}} \cdot {l_{\textrm{pixel}}}$$
where ${n_{\textrm{img}}}$ is the number of pixels of the image in one lateral dimension (assuming a square image), and ${l_{\textrm{pixel}}}$ is the effective camera pixel size (the physical length of one pixel occupied in the specimen). As a result, increasing $NA$ increases the radius of the OTF boundary causing fewer number of constraints imposed on parameters. As cCRLB depends on the number of constraints and the number of pixels within the image, larger $NA$ results in worse cCRLB (Fig. 4). An extreme case will be when the $NA$ is sufficiently large such that OTF occupies the entire field of the Fourier transform of the image. In such a case, there is no constraint imposed on the parameters and cCRLB will be equal to CRLB.

 figure: Fig. 4.

Fig. 4. The relationship between $NA$ and cCRLB. Simulation was conducted on images with 64 by 64 pixels, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. The red circles indicate the OTF boundaries. As $NA$ increases from 0.8 to 1.4 (from left to right), the cCRLB of the image increases.

Download Full Size | PDF

3.4. Pixel size influence on the denoising variance lower bound

Next, we investigated the influence of effective pixel size on the denoising variance lower bound, cCRLB. In order to maintain the size of the field of view as well as the photon flux emitted per unit area, we changed pixel size while adjusting number of pixels of the image and the expected photon count per pixel accordingly. Here, we calculated CRLB and cCRLB in 3 settings: one with an image of 128 by 128 pixels and a pixel size of 40 nm, one with an image of 64 by 64 pixels and a pixel size of 80 nm, and one with an image of 32 by 32 pixels and a pixel size of 160 nm. Since images with smaller pixel sizes result in smaller CRLBs and cCRLBs due to the lower expected photon count per pixel compared to those with larger pixel sizes, to facilitate comparisons among these lower bound maps, we binned both CRLB and cCRLB maps with the size of 128×128 pixels and 64×64 pixels into 32×32 maps (Fig. 5). As shown in Fig. 5, we found that the denoising variance lower bound on image with a pixel size of 40 nm was the lowest indicating the highest level of denoising benefit. Meanwhile, this lower bound was highest when the pixel size was 160 nm suggesting little to no benefit when using denoising algorithms. In contrast, CRLB maps remained equal in all three conditions. This observation shows that choosing a smaller effective pixel size during microscopy experiments, which effectively oversamples a diffraction limited image and therefore increases the number of constraints shared per pixel, benefits the denoising performance with a lower achievable variance bound. The result of NCS denoising algorithm ($\alpha = 1$) performance accorded with our observation (Fig. 5).

 figure: Fig. 5.

Fig. 5. Pixel size influences on denoising variance lower bound and performance of NCS. Specifically, in the case of 160 nm pixel size (left column), there are only 25 pixels outside OTF region forming 25 constraint equations compared to a total of 1024 parameters. As a result, cCRLB is very close to CRLB. In these simulations, the underlying structure, the field of view, and photon count per area were kept the same. Only shot noise was considered. Simulation was conducted with an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and peak expected photon count of 10 per 80 nm by 80 nm pixel. cCRLB maps, CRLB maps, variance and bias maps of NCS denoising algorithm were each calculated from simulated images of 128×128 pixels with a pixel size of 40 nm (binned to 32×32 for comparison), simulated images of 64×64 pixels with a pixel size of 80 nm (binned to 32×32 for comparison) and simulated images of 32×32 pixels with a pixel size of 160 nm, respectively.

Download Full Size | PDF

3.5. Image denoising variance and bias tradeoff

cCRLB represents the theoretical variance lower bound for unbiased estimation of microscopy images. However, to approach unbiased estimation, most denoising algorithms need sacrifice their abilities of reducing estimation variance (Fig. 6). For example, NCS algorithm has an inherent parameter $\alpha $ serving as a balancing parameter between its likelihood function and the prior knowledge. By tuning this parameter, one can balance the tradeoff between bias and variance. We demonstrated this effect in Fig. 6. When $\alpha = 0.1$, we observed the denoised image with relative small bias while its variance performance was far from the theoretical bound (cCRLB). However, as parameter $\alpha $ increased to 3, the variance of denoised image had reached to 7.04 (e-)2 compared to variance of original noisy image 14.89 (e-)2 while introducing the bias with ∼0.28 e- on average at structure regions and the ratio between variance and cCRLB was 1.54 on average. As this parameter further increased ($\alpha = 10$ and $\alpha = 100$), the variance was reduced to 3.55 (e-)2 and 1.70 (e-)2 on average respectively and, at the same time, the absolute bias at the structure regions increased to an average amplitude of 0.35 e- and 0.5 e-. We also included violin plots of variance and bias distributions over a larger range of $\alpha $.

 figure: Fig. 6.

Fig. 6. Trade-off between variance and bias of NCS denoising algorithm. (a) The ideal image and its corresponding cCRLB map. (b) The variance and bias maps calculated using different parameter $\alpha $ in NCS. Every pixel’s cCRLBs and their variances after denoising are shown in a scatter plot (fourth row) with the solid red line representing the cases when cCRLBs and the estimation variances are equal. (c, d) Violin plots of variance and bias distributions on structured regions (regions above average intensity) at different $\alpha $ values. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. Sensor gain map and readout noise ${\sigma ^2}$ map are shown in Fig. 9.

Download Full Size | PDF

4. Conclusion

Our development of cCRLB provides a theoretical variance lower bound for unbiased microscopy denoising algorithms. We demonstrated cCRLB calculation by considering the finite spatial frequency response in a microscope system, representing the common and minimum knowledge available for microscopy imaging denoising. The mathematical framework enables the incorporation prior knowledge either from the imaging system or from the specimen as a set of constraint equations, providing benchmarks for evaluating denoising algorithms. In addition, by using cCRLB, we showed that the readout statistics (readout-noise SD) of individual pixel has limited influence on how precise the pixel intensity can be estimated, suggesting the importance of developing denoising algorithms especially for camera sensors with pixel-dependent readout-noise variations (e.g. CMOS sensor).

Since cCRLB is the variance lower bound for unbiased estimator, for a biased algorithm, it is possible for its performance variance achieving a value lower than cCRLB. In addition, the OTF prior demonstrated in this work represents the common prior knowledge for microscopy image denoising, algorithms that take advantage of other knowledge or assumptions may achieve lower variance than the demonstrated OTF-based cCRLB here. As a result, a new set of constraints must be incorporated when a different type of prior is the major focus of a denoising algorithm such as assumptions on specimen structure [7]. This requires further development. Despite these limitations, we expect our developed method provide a general framework in predicting the precision limit for a proposed microscopy denoising algorithm.

Appendix: Figs. 7 to 10

 figure: Fig. 7.

Fig. 7. Comparison of N2V denoising algorithms trained with different training sets. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy images used as training sets in different structural types. Siemens stars training set contained Siemens stars with 5 to 20 number of spokes, each containing 1000 frames of noisy images. Random block pattern training set was created from 100 different random patterns (resembling blocks of intensities) each containing 1000 corresponding noisy images. Microtubule training set contained 100 simulated different microtubule structures each contacting 1000 corresponding noisy images using worm-like-chain model as described in [3]. Second row shows examples of calibrated noisy image and denoised images from Noise2Void trained with different structure types. Variance map (third row) and bias map (fourth row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps (fifth row) were calculated from variance maps divided by cCRLB pixel by pixel.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. Comparison of three different denoising algorithms at high photon count case. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy image and denoised images from NCS, N2V and ACsN algorithms. Variance map (second row) and bias map (third row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps (fourth row) were calculated from variance maps divided by cCRLB pixel by pixel.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Demonstration of simulation settings and corresponding CRLB and cCRLB. Gain map and variance map were taken from an sCMOS camera (Orca-Flash4.0v3, Hamamatsu Photonics, Japan).

Download Full Size | PDF

 figure: Fig. 10.

Fig. 10. Visualization of matrix ${\boldsymbol P}({\boldsymbol \mu } ){{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$. Simulation was conducted on an image with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. The projected inverse information matrix [Eq. (11)] is 4096 by 4096. The distance between two adjacent strips (shown in zoom-in area) is 64 pixels indicating these two pixels are neighbors vertically. Due to the independent pixel readout, the conventional inverse of information matrix [Eq. (9)] is diagonal. However, after considering the frequency constraint imposed on pixels, the projected inverse of information matrix [Eq. (11)] is no longer diagonal, indicating correlations among pixel estimations during denoising. In another word, the information of one pixel’s expected photon count resides not only in the observation of that pixel but also in the observations of surrounding pixels.

Download Full Size | PDF

Funding

National Institute of General Medical Sciences (R35GM119785).

Acknowledgments

We would like to thank Fan Xu for suggestions on our simulations, fruitful discussions and helps on revising the manuscript. We would also like to thank Benjamin Brenner for suggestions on our simulations and manuscript.

Disclosures

The authors declare no conflicts of interest.

References

1. B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (john Wiley & sons, 2019).

2. A. Krull, T. O. Buchholz, and F. Jug, “Noise2void-Learning denoising from single noisy images,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR) (2019), Vol. 2019-June, pp. 2124–2132.

3. S. Liu, M. J. Mlodzianoski, Z. Hu, Y. Ren, K. McElmurry, D. M. Suter, and F. Huang, “sCMOS noise-correction algorithm for microscopy images,” Nat. Methods 14(8), 760–761 (2017). [CrossRef]  

4. J. Boulanger, C. Kervrann, P. Bouthemy, P. Elbau, J. B. Sibarita, and J. Salamero, “Patch-based nonlocal functional for denoising fluorescence microscopy image sequences,” IEEE Trans. Med. Imaging 29(2), 442–454 (2010). [CrossRef]  

5. B. Mandracchia, X. Hua, C. Guo, J. Son, T. Urner, and S. Jia, “Fast and accurate sCMOS noise correction for fluorescence microscopy,” Nat. Commun. 11(1), 94 (2020). [CrossRef]  

6. A. Buades, B. Coll, and J. M. Morel, “A non-local algorithm for image denoising,” in Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR) (2005), Vol. II, pp. 60–65.

7. M. Weigert, U. Schmidt, T. Boothe, A. Müller, A. Dibrov, A. Jain, B. Wilhelm, D. Schmidt, C. Broaddus, S. Culley, M. Rocha-Martins, F. Segovia-Miranda, C. Norden, R. Henriques, M. Zerial, M. Solimena, J. Rink, P. Tomancak, L. Royer, F. Jug, and E. W. Myers, “Content-aware image restoration: pushing the limits of fluorescence microscopy,” Nat. Methods 15(12), 1090–1097 (2018). [CrossRef]  

8. A. Levin and B. Nadler, “Natural image denoising: Optimality and inherent bounds,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR) (2011), pp. 2833–2840.

9. J. W. Goodman, Introduction to Fourier Optics 3ed (Roberts and Company Publishers, 2005).

10. K. R. Rao and P. C. Yip, The Transform and Data Compression Handbook (CRC press, 2018).

11. F. Huang, T. M. P. Hartwich, F. E. Rivera-Molina, Y. Lin, W. C. Duim, J. J. Long, P. D. Uchil, J. R. Myers, M. A. Baird, W. Mothes, M. W. Davidson, D. Toomre, and J. Bewersdorf, “Video-rate nanoscopy using sCMOS camera-specific single-molecule localization algorithms,” Nat. Methods 10(7), 653–658 (2013). [CrossRef]  

12. R. J. Ober, S. Ram, and E. S. Ward, “Localization Accuracy in Single-Molecule Microscopy,” Biophys. J. 86(2), 1185–1200 (2004). [CrossRef]  

13. S. M. Kay, Fundamentals of Statistical Signal Processing (Prentice Hall PTR, 2013).

14. J. D. Gorman and A. O. Hero, “Lower Bounds For Parametric Estimation with Constraints,” IEEE Trans. Inf. Theory 36(6), 1285–1301 (1990). [CrossRef]  

15. H. Cramer, Mathematical Methods of Statistics. (Princeton University Press, 1999).

16. C. R. Rao, “Information and the Accuracy Attainable in the Estimation of Statistical Parameters,” in Breakthroughs in Statistics (Springer, 1992), pp. 235–247.

17. K. B. Petersen and M. S. Pedersen, The Matrix Cookbook (Technical University of Denmark, 2008).

18. Y. Li, S. Liu, and F. Huang, “cCRLB source code,” github.com (2020), https://github.com/HuanglabPurdue/cCRLB.

Supplementary Material (1)

NameDescription
Code 1       Source code to calculate cCRLB

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1.
Fig. 1. Schematic diagram of obtaining cCRLB, bias, and variance maps for evaluating performance of microscopy image denoising algorithms. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. CRLB is the variance lower bound on estimation of ideal image calculated from the noise model [Eq. (4)–(10)]. OTF is the optical transfer function. cCRLB is the variance lower bound on estimation of ideal image while taking CRLB and additional prior knowledge on frequency constraint [Eq. (3)] into consideration. Bias map and variance map (on the right) were calculated pixel by pixel from 100 frames of denoised images (using NCS denoising algorithm).
Fig. 2.
Fig. 2. Comparison of three denoising algorithms and cCRLB. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy image (raw noisy image subtracted by the offset and divided by the gain map) and denoised images from NCS ( $\alpha = 3$ ), N2V and ACsN algorithms. Variance map (second row) and bias map (third row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps are shown in the fourth row.
Fig. 3.
Fig. 3. Comparison of CRLB and cCRLB at different readout noise levels. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm (a) The ideal image with the target pixels marked. (b) Comparisons between CRLB and cCRLB with respect to increasing readout-noise SD in the target pixel while those of the remaining pixels were kept constant. (c) The ratios of CRLB/cCRLB of the target pixel and the mean ratios of its neighbors (8 surrounding pixels) at different readout noise levels of the target pixel.
Fig. 4.
Fig. 4. The relationship between $NA$ and cCRLB. Simulation was conducted on images with 64 by 64 pixels, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. The red circles indicate the OTF boundaries. As $NA$ increases from 0.8 to 1.4 (from left to right), the cCRLB of the image increases.
Fig. 5.
Fig. 5. Pixel size influences on denoising variance lower bound and performance of NCS. Specifically, in the case of 160 nm pixel size (left column), there are only 25 pixels outside OTF region forming 25 constraint equations compared to a total of 1024 parameters. As a result, cCRLB is very close to CRLB. In these simulations, the underlying structure, the field of view, and photon count per area were kept the same. Only shot noise was considered. Simulation was conducted with an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and peak expected photon count of 10 per 80 nm by 80 nm pixel. cCRLB maps, CRLB maps, variance and bias maps of NCS denoising algorithm were each calculated from simulated images of 128×128 pixels with a pixel size of 40 nm (binned to 32×32 for comparison), simulated images of 64×64 pixels with a pixel size of 80 nm (binned to 32×32 for comparison) and simulated images of 32×32 pixels with a pixel size of 160 nm, respectively.
Fig. 6.
Fig. 6. Trade-off between variance and bias of NCS denoising algorithm. (a) The ideal image and its corresponding cCRLB map. (b) The variance and bias maps calculated using different parameter $\alpha $ in NCS. Every pixel’s cCRLBs and their variances after denoising are shown in a scatter plot (fourth row) with the solid red line representing the cases when cCRLBs and the estimation variances are equal. (c, d) Violin plots of variance and bias distributions on structured regions (regions above average intensity) at different $\alpha $ values. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. Sensor gain map and readout noise ${\sigma ^2}$ map are shown in Fig. 9.
Fig. 7.
Fig. 7. Comparison of N2V denoising algorithms trained with different training sets. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy images used as training sets in different structural types. Siemens stars training set contained Siemens stars with 5 to 20 number of spokes, each containing 1000 frames of noisy images. Random block pattern training set was created from 100 different random patterns (resembling blocks of intensities) each containing 1000 corresponding noisy images. Microtubule training set contained 100 simulated different microtubule structures each contacting 1000 corresponding noisy images using worm-like-chain model as described in [3]. Second row shows examples of calibrated noisy image and denoised images from Noise2Void trained with different structure types. Variance map (third row) and bias map (fourth row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps (fifth row) were calculated from variance maps divided by cCRLB pixel by pixel.
Fig. 8.
Fig. 8. Comparison of three different denoising algorithms at high photon count case. Simulation was conducted on images with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. First row shows examples of calibrated noisy image and denoised images from NCS, N2V and ACsN algorithms. Variance map (second row) and bias map (third row) were calculated using 100 frames of corresponding noisy or denoised images. Variance/cCRLB ratio maps (fourth row) were calculated from variance maps divided by cCRLB pixel by pixel.
Fig. 9.
Fig. 9. Demonstration of simulation settings and corresponding CRLB and cCRLB. Gain map and variance map were taken from an sCMOS camera (Orca-Flash4.0v3, Hamamatsu Photonics, Japan).
Fig. 10.
Fig. 10. Visualization of matrix ${\boldsymbol P}({\boldsymbol \mu } ){{\boldsymbol I}^{ - 1}}({\boldsymbol \mu } )$ . Simulation was conducted on an image with 64 by 64 pixels, an objective $NA$ of 1.4, an emission wavelength in vacuum of 700 nm and a pixel size of 80 nm. The projected inverse information matrix [Eq. (11)] is 4096 by 4096. The distance between two adjacent strips (shown in zoom-in area) is 64 pixels indicating these two pixels are neighbors vertically. Due to the independent pixel readout, the conventional inverse of information matrix [Eq. (9)] is diagonal. However, after considering the frequency constraint imposed on pixels, the projected inverse of information matrix [Eq. (11)] is no longer diagonal, indicating correlations among pixel estimations during denoising. In another word, the information of one pixel’s expected photon count resides not only in the observation of that pixel but also in the observations of surrounding pixels.

Equations (20)

Equations on this page are rendered with MathJax. Learn more.

[ F T { U } ] p , q = j = 0 n 1 k = 0 n 1 U j , k e i 2 π n ( p j + q k ) = j = 0 n 1 k = 0 n 1 W p , j U j , k W k , q = [ W U W ] p , q = 1 p T W U W 1 q
W p , q = e i 2 π p q n .
E p , q = ( 1 p T W U W 1 q ) T ( 1 p T W U W 1 q ) = 1 q T W T U T W T 1 p 1 p T W U W 1 q = 0
μ = [ μ 1 , μ 2 , , μ N ] T
z = [ z 1 , z 2 , , z N ] T
Z j = g j S j + Y j + o j
p Z j ( z j | μ j ) = τ = 0 p S j ( τ | μ j ) p Y j ( z j g j τ o j ) = τ = 0 e μ j μ j τ τ ! 1 2 π σ j exp { ( z j g j τ o j ) 2 2 σ j 2 }
Σ ( μ ) = E [ ( μ ^ μ ) T ( μ ^ μ ) ]
Σ ( μ ) I 1 ( μ )
[ I ( μ ) ] j , j = 1 μ j 2 [ τ = 0 τ p S j ( τ μ j ) p Y j ( z j g j τ o j ) ] 2 τ = 0 p S j ( τ μ j ) p Y j ( z j g j τ o j ) d z j 1
Σ ( μ ) P ( μ ) I 1 ( μ )
P ( μ ) = I I 1 ( μ ) G T ( μ ) [ G ( μ ) I 1 ( μ ) G T ( μ ) ] 1 G ( μ ) .
G ( μ ) = [ G 1 ( μ ) , G 2 ( μ ) , , G K ( μ ) ] T = 0
[ G ( μ ) ] j , k = G j ( μ ) μ k .
G j ( μ ) = vec ( E p , q U ) T .
[ E p , q U ] j , k = E p , q U j , k .
D = W T 1 p 1 p T W , b = W 1 q , D = b b T = W 1 q 1 q T W T .
E p , q = b T U T D U b
E p , q U = 2 D U b b T = 2 D U D
R OTF = 2 N A λ n img l pixel
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.