1 Introduction
Using fluorescence microscopy with singlemolecule sensitivity, it is now possible to track the movement of individual fluorophoretagged molecules such as proteins and lipids in the cell membrane with nanometer precision. This ability has been used to characterize the diffusional properties of molecules Saxton97:SingleParticleTracking, as well as monitor their organization relative to other molecules Dunne09:SMF and cellular structures Andrews08:SingleParticleTracking, with the potential to work below the diffraction limit of light.
The design of automatic, efficient computer algorithms for data analysis is an extremely important facet of such experiments. Ideally, the program should be able to tolerate the low signal to noise ratios inherent to singlemolecule data, must work at high object density and be able to cope with large data volumes at reasonable speed. Suitable methods can usually be decomposed into several sequential steps: filtering of the data, identification of fluorescent objects and their precise positions, and ‘tracking’ of these objects by linking together their positions over an image sequence. For precise tracking care must be taken to optimise each of these stages.
Previous work in Single Molecule Fluorescence image processing has focussed on the tracking step of the problem, with several suggested approaches resulting in excellent performance, including a grid algorithm Bruckbaur07:SMF, and several Monte Carlobased algorithms Yoon08:SMF; Dunne09:SMF. However, with one exception Yoon08:SMF, the initial denoising step utilized conventional denoising algorithms, such as Gaussian filtering. Afterwards, the detection algorithms are applied to the denoised images in order to separate signals and background noise. As a result, the tracking algorithms had to assume a cluttered environment with false alarms and missed detections, impacting their performance. Thus in this paper we concentrate on the first step in the data analysis procedure: ‘data filtering’, or ‘denoising’.
Denoising is commonly carried out using one of two approaches: Gaussian smoothing or wavelet denoising. In the first, the image is smoothed through convolution with a Gaussian or average kernel, which is related to matched filtering for additive uncorrelated noise Turin60:MatchedFilter. Use of a Gaussian filter is justified since the fluorescence profile of an individual molecule is well approximated by (and hence frequently fitted to) a Gaussian function Zhang07:GA4Fluorescence. In cases where the noise addition is linear or stationary, a Wiener filter is also used to reduce the noise of the image Wiener64:WienerFilter. For the nonstationary problem, wavelet based denoising algorithms are also wellknown for image restoration Donoho93:nonlinearwavelet; Marin02:SpotDetection.
In this work we develop and describe a new denoising algorithm based on a Gaussian Markov Random Field (GMRF) model. Our proposed algorithm is a fully Bayesian approach with few tuning parameters: we need to set only a small relative threshold parameter and specify hyperparameters to build prior distributions. The performance of this method relative to previous methods is evaluated using both synthetic and real single molecule data and found to display significant advantages.
The paper is structured as follows. We first introduce the Gaussian Markov random field in a latent Gaussian model (Section 2) which is then used to propose mathematical models for denoising images (Sections 3 and 4). Finally, we compare the results of the various algorithms when applied to synthetic and real data (Section 5).
2 Intrinsic Gaussian Models
Gaussian Markov random fields (GMRFs) are Gaussian fields defined on a discrete grid with a Markov property of conditional independence of a component with all others given its neighbours Mardia88:GMRF; GMRFbook:Rue. They have seen widespread application in statistical modelling, for example in spatiotemporal models Besag91:GMRF and dynamic linear models West97:LDM.
In this paper, a GMRF
is an one dimensional vector which corrsponds to a two dimensional image lattice
GMRFbook:Rue. Let be differenced values of at a site on the lattice with neighbours. There are many possible ways to define a GMRF for through the . For example, a first order model can be defined through assuming the differences with the 4 nearest neighbours ,to be Gaussian with mean 0 and precision (inverse of variance)
. The distribution of is of multivariate Gaussian form:(1)  
(2)  
(3) 
where the precision matrix (inverse of covariance matrix) , and is a matrix with elements , if the pixels represented by the th and th components are neighbours, and 0 otherwise. The resulting precision matrix
is not of full rank, in which case it is known as an intrinsic GMRF (IGMRF). This form is used often as a prior in Bayesian inference. For denoising, it has the desirable properties of placing more prior weight on smooth images while avoiding the specification of a mean value of
. Under very general conditions, the posterior of will become a proper distribution GMRFbook:Rue.3 The Denoising Model
Note that we focus on denoising a single image in this paper although we have a sequence of images for processing. For simplicity, we assume images in the sequence are independent so that the denoising algorithm can be applied separately. The dependency of spots across the sequence of images is accounted for in the tracking process.
3.1 Linear Model for an intrinsic GMRF
We assume a Gaussian model for an observed image in terms of the underlying signal and some regression coefficients :
(4) 
where for an unstructured term and
denotes the Gaussian distribution. In this model,
are lattice location indices so that models any global background noise trend. Therefore, the set of parameters are denoted .3.2 Priors
We assume conjugate prior forms of the regression coefficients
(5) 
where , giving a weak prior, and
(6)  
(7) 
where
denotes the gamma distribution with hyperparameters
and .According to the linearity of the hidden parameters in Eq. (4), we can apply RaoBlackwellization technique Casella96:RaoBlackwellisation to reduce the dimension of the hidden variables and variances. Now, we have marginalized likelihood defined by
(8)  
(9)  
(10) 
where .
3.3 Problems with the IGMRF prior
GMRFbook:Rue describes several limitations of the intrinsic GMRF. The limitations will bring following problems in denoising images:

Unwanted blurring effects can be introduced because the IGMRF penalises discrete boundaries or sharp gradients between neighbours;

Weak signals can be ignored by smoothing with the background;

Stationary (homogeneous) fields across an image rarely exist in practice.
In spite of its convenient mathematical properties, these three concerns present difficulties for this application. A variant of the conventional IGMRF is therefore introduced that assigns different interaction weights to pixel neighbours according to whether the pixel is classified as being part of a spot (the molecule of interest is present) or the background. We name this variant of the IGMRF as the Heterogeneous IGMRF (HIGMRF).
3.4 Heterogeneous IGMRF (HIGMRF)
In this application, pixels in the image can be classified as either a spot or background. This provides a way to allow for nonstationarity in an IGMRF. For pixel let if it is a spot and if it is background. The conditional distribution is of the same Gaussian form of Eq. (3), and the precision matrix is defined through differences as with the IGMRF, but the differences are weighted by between a background pixel and neighbouring spot pixels, so
(11) 
for
. The effect of this prior is to give higher probability to images with large contiguous background regions whilst preventing oversmoothing of spots. The value of
must be specified; is used here and works well for the images used in this paper. This model is related to the Ising model commonly used in image processing Gelman84:GibbsSampling. (Although thecan be estimated from the
, we do not estimate in Gibbs scheme because it is time consuming.4 Posterior Calculations
The unknown parameters in the marginalized form are and and so the goal is to evaluate . This is implemented by a Gibbs sampler that samples from
(12) 
where and
Then and are sampled separately from . However, it is not straightforward to draw samples from the marginalized conditional posterior since there is no exact conditional distribution which we can easily generate samples. Thus, we draw the from the complete posterior (which is a sort of Data Augmentation technique) and then we generate samples on the model parameters and .
(13) 
where and .
Now, we can generate samples of and separately from , which are of the form:
(14)  
(15) 
where
(16)  
(17)  
(18)  
(19) 
4.1 Sampling in the HIGMRF case
When is an HIGMRF, it depends on . Good performance of the algorithm can be obtained if is updated at each iteration of the Gibb’s sampler by a locally defined threshold ; if (background) and otherwise (spot), reflecting the fact that the background is darker. The details of the algorithm are described in algorithm 1. The threshold is specified through a parameter ; empirically it has been found that works well for the images used in this paper. More details to find such a hard threshold will be referred to in discussion section.
4.2 Algorithm for HIGMRF
The Gibb’s sampler is initialized with , defined as in the usual IGMRF case of Section 2 . The th iteration of the Gibb’s sampler is:

is sampled from .

is sampled from .

is sampled from the marginalized posterior .

If the HIGMRF prior is used for :

is calculated from from algorithm 1.

is calculated from .

The denoised image is reported as the sample average of the , following an initial burnin period.
(a) Ground  (b) Gaussian  (c) Average  (d) Wiener  (e) Wavelet  (f) LGW 
Truth  (Ga)  (Av)  (Wi)  (Wav)  
(g) Noisy image  (h) ROF  (i) NLM  (j) TV  (g) IGMRF  (h) HIGMRF 
4.3 Convergence Diagnostics
Convergence of the MCMC method is checked through the potential scale reduction factor (PSRF) diagnostic of Gelman92:MCMCconvergence that monitors selected outputs of the MCMC algorithm and is perhaps more accurately described as output analysis Brooks98:MCMCconvergence; Brooks98:MCMCconvergence2. For scalar valued sequences of length , , , we define the following quantities:
(20)  
(21) 
where . The is defined to be
The is computed for all univariate parameters separately following burnin. Convergence is consistent with values below to 1.2 Gelman92:MCMCconvergence. If this criterion is not satisfied for all parameters then it is indicative of lack of convergence and the simulation run needs to be lengthened.
4.4 Evaluation of the denoising performance
Let and be the estimated image by a particular denoising algorithm and its ground truth image respectively.
4.4.1 Root Mean Square Error (RMSE)
The most common and simplest metric for the image quality assessment is Root Mean Square Error (RMSE) which is computed by averaging the squared intensity differences of and by .
4.4.2 Peak SignaltoNoise Ratio (PSNR)
Peak SignaltoNoise Ratio (PSNR) is a wellknown quality metric which is an extension of the RMSE by where is the maximum intensity of the estimated image pixels.
4.4.3 Kullback Leibler Distance (KLD)
Let and be the true distribution and the estimated distribution respectively by where the distributions are estimated via quantification into bins. In this paper, we set .
4.4.4 Structural SIMilarity (SSIM)
The other metric is structural similarity based image quality assessment which is able to explain the strong dependencies between pixels of the highly structured images Wang04imagequality.
(22) 
where , , , , and . In this paper, for simplicity we use which is commonly defined as ”universal quality index (UQI)” although this produces unstable results when the denominator of Eq. (22) is very close to zero Wang02:UQI.
5 Application to synthetic and real single molecule data
5.1 Synthetic images
Initially we sought to quantitatively compare the performance of the various algorithms using synthetic data. Synthetic data consisted of 50 randomly generated groundtruth images of a pixel region with a small number of single molecule signals, to which varying levels of noise were added (resulting in signal to noise ratios in the range of 510 dB for different noise levels). The various parameters used in generation of this data are given in Table I. Next, various filtering algorithms were applied to the data. Figure 1 shows an example of simulated data (Fig. 1 a and g) and the denoised images obtained through different algorithms: Gaussian filter (Ga), Average filter (Av), Wiener filter (Wi), Wavelet filter (Wav) Donoho93:nonlinearwavelet, Log Garbor Wavelet (LGW) based filter Kovesi_phasepreserving, Total Variation approach by Rubin et al. (ROF) Rudin92:TV, NonLocal Mean filter (NLM) Buades05:anonlocal, Total Variation by Split Bregman Method (TV) Goldstein09:Denoising, IGMRF, and HIGMRF of Fig. 1. By eye it appeared that the HIGMRF filtered images came closest to matching the underlying groundtruth. Analyzing the root mean square error (RMSE), Kullback Leibler Distance (KLD), Peak SignaltoNoise Ratio (PSNR) and Structural SIMilarity (SSIM) between the intensities of the ground truth relative to recovered intensities in the filtered images (Fig. 2) confirmed this supposition, showing that intensities recovered after HIGMRF were much closer to the truth than the other methods.
parameter  value  parameter  value 

(a) RMSE  (b) KLD 
(c) PSNR  (d) SSIM 
5.2 Single molecule fluorescence images
The performance of the algorithms was then tested with real single molecule fluorescence data, obtained on living human Tcells, a cell of the immune response. To create a suitable sample, the immunoglobulin protein CD86 was genetically introduced into human Tcells so that it was present at the cell membrane at low concentrations, as described in Dunne et al.Dunne09:SMF. These molecules were labelled with an organic dyelabelled antibody fragment specific to CD86 (Alexa Fluor 488 labelled Fab molecule, derived from the antiCD86 antibody BU63, see James et al.James07:SingleMolecule for details), and imaged using Total Internal Reflection Fluorescence Microscopy, as described previously Dunne09:SMF.
For this real experimental image, the performances of ten different denoising algorithms as shown in synthetic dataset were compared for one frame (Fig. 3 a). The results from these filtering steps are shown in Fig. 3 (b)(k). The figures show only chopped ared ( size) to clearly see the denoising effect. Compared to other approaches the HIGMRF preserves the spot signals of interest while it suppresses the global and local noises effectively. The convergence of our Gibbs sampler during IGMRF and HIGMRF filtering of this data was checked using the PSRF diagnostic method described in Section 4.3. The PSRF was found to converge to a value below 1.2 for all components of (Fig. 4).
(a) Observation  (b) Gaussian  (c) Average  (d) Wiener  (e) Wavelet  (f) LGW 
(Ga)  (Av)  (Wi)  (Wav)  
(g) ROF  (h) NLM  (i) TV  (j) IGMRF  (k) HIGMRF 
(a) Histograms of sampled values of  (b) The convergence of 
6 Discussion
Several issues for the current proposed algorithm remain to be resolved and discussed. First of all, inferring is a key step for the proposed algorithm. Basically, is regarded as a set of clustered labels given continuous values . Therefore, we can apply any clustering algorithms with two clusters. For instance, we can partition (cluster) into different clusters by using means algorithms to infer with . However, it is rather difficult to build locality based clusters and it might infringe the properties of Gibbs sampling. The best way to reconstruct is to generate them in a Gibbs sampling scheme (that is, singlespin updates which is called ’heat method’ in the statistical physics). However, it takes a lot of time to reach the convergence since the dimension of is too large to obtain sufficiently good results in a practical time. (As we know, Gibbs sampling requires a myriad of samples to restore joint posterior when variables (dimensions) are highly dependent.) Worse, the variables of
are not continuous variables but binary variables, which are often rather difficult to infer. We can apply the SwendsenWang model which is known as relatively fast sampling algorithm for Ising model by adding auxiliary variable for bonds. However, inferring
via Monte carlo sampling even with such a fast SwendsenWang model is still problematic since it brings a lot of time consuming tasks for the convergence. Therefore, we adopt a simple threshold algorithm in order to reduce time complexity although the threshold based algorithm has the sensitivity problem. The threshold based algorithm can be regarded as a sort of a simple clustering algorithm. Also, it works as if the proposal function is a Dirac delta function in Gibbs sampling. From this point of view, the simple threshold approach is practically robust and useful. In this current version, we set a very low value by in order to allow relatively dim spot regions to be classified. The threshold may be estimated in the computation as well, although the fixed value worked well in the examples described here. As a further work, one might consider assuming unknown and try to infer its value in the Gibb’s sampler scheme. In addition, we may be able to jointly generate samples of andon by embedding Bernoulli distribution into the system. This joint estimation will improve our current algorithm a lot since it is known that such a highly dependent parameters can be improved by blocked Gibbs sampling technique.
In addition, although most of the hyperparameters in the model — and for — can be assigned vague priors, an informative prior on and is necessary. These are the hyperparameters for the precision of the (H)IGMRF and the solution is sensitive to their value.
7 Conclusion
In this paper, a hierarchical model has been proposed to denoise single molecule fluorescence images, in which a Gaussian Markov random field (GMRF) is used as a well defined prior of a latent Gaussian model. A variant of the intrinsic GMRF, called the heterogeneous intrinsic GMRF (HIGMRF), was introduced in order to maintain the attractive properties and computational advantages of the intrinsic GMRF while addressing some of its weaknesses, such as blurring and stationary. These are obtained by assigning nonstationary interaction weights to pairs of pixels of the image. The weights and hence precision matrix (inverse of covariance matrix) for the HIGMRF are updated sequentially in a Gibbs framework. It has been demonstrated that the use of the HIGMRF reduces global and local noise from images effectively.
Comments
There are no comments yet.