Directional bilateral filters for smoothing fluorescence microscopy images

In biological imaging, the working distance of the objective lens determines the depth to which one can image into the sample.⁴¹ It is also well known that numerical aperture (NA) of the objective lens decreases as the working distance increases and the resolution of an optical microscope is determined by the NA.^20,38 As a result of this trade off, the images obtained using low-NA objectives have less contrast and poor resolution. In order to address these problems, a number of clearing agents have been developed.¹⁷ By employing clearing agents, it becomes possible to identify the structure at a cellular resolution but live imaging of the biological specimens is not possible. Moreover, such techniques require customized objectives for clearing agents.^13,17 In this context, post-processing of the images acquired with a low-NA objective lens plays an important role.

From Fig. 1, we observe that images obtained at low NA are noisy and the filaments do not appear continuous. In the presence of noise, it becomes difficult to resolve two or more filaments that are close to each other and hence there is a need to smooth without significantly affecting the resolution. The filaments contain important structural information and are directional. It is crucial to retain and emphasize the directionality information while suppressing noise.

In conventional widefield fluorescence imaging, the entire sample is illuminated at a specific wavelength. Every fluorophore in the specimen acts as an independent source of light and contributes to the image formed at the detector. Along with the photons captured from the focal plane, the detector also collects the fluorescence photons from the planes that are above and below the focal plane. This results in image degradation and blurring of the in-focal details. These problems were solved with the advent of confocal fluorescence microscopy.^26,28 It employs a pinhole in front of the detector to eliminate the out-of-focus photons and allows only the in-focus photons to impinge on the detector. This makes the confocal microscope a serial acquisition system than a parallel one and it requires a bright source of light for illumination.^16,29 With a smaller pinhole diameter, it is possible to get sharper images but the photon count goes down, which reduces the signal-to-noise ratio (SNR) in the acquired images.³² On the other hand, a pinhole with a larger diameter introduces additional statistical noise originating from the fluorophores outside the plane of interest, autofluorescence of the specimen and also from the optical components.⁹ The integrated effect limits the image contrast.

In this paper, we deal with noise suppression in fluorescence images that contain important directional information. More specifically, we have obtained low NA images of F-actin filaments using a widefield epifluorescence microscope and a confocal microscope. F-actin is the protein polymer responsible for the formation of cytoskeleton. It is bountiful in all Eukaryotic cells. The richest area of actin filaments in a cell lies in a narrow region just below the plasma membrane. The physical diameter of these filaments is approximately 7nm. F-actin filaments are widely studied in cell biology because of their essential role in numerous cellular processes like phagocytosis, morphogenesis, cytokinesis, and endocytosis.³⁹ Dynamic behavior of the cytoskeleton actin also governs cell shape and cell motility.³⁹ Without cytoskeleton, cells will not be able to resist physical stress. It is also widely studied in pathology, because the abnormalities in actin dynamics play a key role in cancer, heart diseases, vascular diseases, myofibrillar myopathies, neurological disorders,³⁰ etc.

This paper is organized as follows: In Section II, we compare linear filtering with bilateral filtering. We also review some related work on bilateral filters, discuss the importance of selecting the optimal parameters, and the methodology behind choosing the parameters. In Section III, we propose a modification to the bilateral filters to enhance edge preservation by taking into account the local orientation. The optimal selection of the filter parameters is discussed in Section IV. Results obtained using the proposed technique on images acquired through fluorescence microscopy and comparisons are shown in Section V.

Linear filtering employs a space-invariant, linear kernel such as a Gaussian to perform uniform smoothing and suppress noise. To evaluate the filtered output at a pixel of interest p, which is the ordered pair of x and y coordinates, the averaging is localized to an Ω-neighborhood of p. Let q be a pixel in this neighborhood and let its intensity be y_q. The output of the filter, given the kernel ϕ_p,q is given by

where h_p is the normalizing factor given by $h p = ∑ q ∈ Ω ϕ p , q$⁠. A Gaussian domain kernel is given by

where ‖p − q‖ is the geometric distance between the pixels p and q. The Gaussian kernel assigns higher weights to pixels that are close to p, and the weights decay, according to the decay parameter σ_d, as the distance increases. Since this kernel acts only on the domain, that is, it is a function of only the distance between the two pixels (the intensities of the pixels are not considered), it is referred to as a domain kernel.

However, high-frequency components such as edges are also smoothed when such a lowpass filter is used. Since edges contain the primal sketch of an image, it is desirable to preserve them. Nonlinear filters that are data-adaptive were designed to smooth images without blurring the edges. Anisotropic diffusion described by Perona and Malik³⁴ was first used to achieve edge-preserving smoothing. Subsequently, Aurich and Weule employed nonlinear modifications of Gaussian filters.³ Tomasi and Manduchi proposed generalized bilateral filters whose range kernel suppresses outliers from smoothing to achieve edge preservation.⁴⁰ The bilateral filter ϕ is obtained by combining a domain kernel and a range kernel

where the domain kernel w_p−q depends on the geometric distance between the pixel of interest p and a neighboring pixel q. It assigns weights that fall off with increasing geometric distance. The range kernel r(y_p − y_q) measures the similarity between the intensity at the pixel of interest y_p and a neighbourhood pixel y_q. It assigns higher weights to pixels that have relatively similar intensities. The bilateral filter output is

where h_p is a normalizing factor given by

Elad showed that the bilateral filter and anisotropic diffusion emerge from a Bayesian framework.¹⁴ 

The Gaussian bilateral filter (GBF) has Gaussian domain and range kernels and is given by

The parameters σ_d and σ_r control the rates at which the Gaussian functions decay and directly determine the degree and locality of smoothing.

The edge-preserving smoothing property of the bilateral filter over the standard Gaussian filter is illustrated in Figures 2 and 3.

Optimizing the parameters of the bilateral filter is crucial in maximizing its performance for image denoising. Consider an image x (vector representation of an image) corrupted by additive white Gaussian noise (AWGN) b of zero-mean and σ²I covariance matrix. The noisy image y is given by y = x + b. The denoised image $x ˆ =f ( y )$ should be an accurate estimate of x. Here, f(y) = [f_n(y)]_1≤n≤N is an N-dimensional vector function and N is the total number of pixels in the image. The mean-squared error (MSE), which quantifies the closeness of the filtered image to the original, is defined as

The ground truth x is not available in practice and hence, it is not feasible to compute the MSE. The ground-truth-dependent terms in (7) are written as

The term $E { f T ( y ) x }$ is evaluated as shown by Blu and Luisier.⁷ For the sake of completeness, we provide the derivation here. For more details, we refer to the original paper by Blu and Luisier.⁷ The key ingredient is a property of the Gaussian probability density q(b_n): b_nq(b_n) = − σ²q′(b_n). Denoting the partial expectation over the nth component of the noise by $E b n { ⋅ }$⁠, we have the following:

Taking the expectation over the remaining components of the noise gives

Since the expectation is a linear operator, we get

That is, the expressions f^T(y)x and f^T(y)y − σ²div{f(y)} have the same expectation. It follows that

is an unbiased estimate of the MSE.

Kishan and Seelamantula²² and Peng and Rao³³ derived SURE for the Gaussian bilateral filter. The differential of the Gaussian filter output with respect to the noisy image (required for the calculation of the divergence term) is obtained as

For the Gaussian domain kernel, the derivative of (6) with respect to y_p is

In arriving at (9), we have assumed prior knowledge of the noise variance σ². However, in practice, it has to be estimated. When images are corrupted by additive white Gaussian noise, it is estimated as²⁴ 

where h is a 2-D highpass filter specified by the mask

The noise is estimated from the highpass subbands. The median effectively suppresses outliers in the highpass subband, which correspond to edges in the image.

Due to the nonlinearity and shift-variance introduced by the range kernel, bilateral filtering is computationally expensive in its standard form. A number of solutions have been proposed to accelerate them. Paris and Durand derived criteria for downsampling in space and intensity to come up with a fast approximation of the bilateral filter.³¹ A constant-time algorithm for fast bilateral filtering has been proposed by Porikli³⁵ and Chaudhury et al.¹⁰ Yang et al. achieved substantial acceleration at the cost of quantization.⁴⁵ 

Modifications of the bilateral filter have found widespread use in a number of image processing tasks such as denoising,^5,11 illumination compensation,¹⁵ optical-flow estimation,⁴⁴ demosaicing,³⁶ edge detection,²¹ etc. Structure-oriented Gaussians have been used for performing detail-preserving smoothing in seismic signals.⁴² Awate and Whitaker developed a technique that constructs a filtering strategy based on the image statistics learned from the degraded data.⁴ 

Considerable work has also been done on optimizing the parameters of the bilateral filter for improving denoising performance. Peng and Rao used Stein’s unbiased risk estimate (SURE) to find the optimal parameters of the Gaussian bilateral filter.^33,37 Kishan and Seelamantula achieved this goal for a bilateral filter with a raised-cosine range kernel.²² Chen and Shu used Chi-square unbiased risk estimate (CURE) for optimizing bilateral filter parameters in squared magnitude MR images.^12,23 Anand and Sahambi used an approach based on power-law scaling of the inverse of local statistics for pixel-wise estimation of the range parameter.² Hashii et al. used the Hellinger distance to measure the similarity between the estimated noise distribution and the assumed noise distribution to tune the filter parameters.¹⁸ 

Traditionally, bilateral filters have employed fixed domain kernels. Our intuition is that if the domain kernel can be locally controlled and adapted to smooth along the dominant local orientations in the image, the influence of outliers can be suppressed while smoothing and directional selectivity can be retained. We design an oriented domain kernel, where the orientation and anisotropy is estimated by means of a structure tensor.⁶ The domain kernel is combined with a range kernel to preserve edges. The optimal parameters of the directional bilateral filter are obtained by minimizing the SURE cost. We evaluate SURE for the directional bilateral filter and show that SURE follows the MSE closely. We show that directional bilateral filter performs orientation-preserving smoothing.

The standard Gaussian domain kernel is symmetric around the center p of the window. Since we want to incorporate orientation and anisotropy of image structures while smoothing, we use an oriented Gaussian domain kernel. The anisotropic domain kernel is given by

where

(m_p, n_p), (m_q, n_q) are the coordinates of the pixels p and q, respectively. Figure 4 shows an example of an anisotropic Gaussian domain kernel.

The additional parameters γ₁, γ₂, and θ control the scaling and the orientation of the Gaussian. They allow for smoothing along a particular direction by taking into consideration the orientation and anisotropy. We compute the parameters γ and θ locally using the structure tensor approach so that the Gaussian domain kernel can be steered accordingly (cf. Figure 5).

A structure tensor gives accurate orientation estimate and anisotropy measure in local neighbourhoods. Let the grayscale image be denoted by I. The difference of Gaussians (DoG) kernel is used to compute the gradient of the image $∇I= I x I y$⁠. The 2-D structure tensor J_ρ is a smoothed version of the second moment matrix (∇I)(∇I)^T. Smoothing is performed by convolving the matrix components with a Gaussian kernel G_ρ, ρ denoting the standard deviation:

J_ρ is a symmetric, positive semidefinite matrix. The information about orientation and anisotropy is obtained by eigenvalue decomposition. The eigenvalues are obtained directly from J_ρ as:⁶ 

The parameters θ and γ are obtained as follows:

1. The domain kernel is oriented along the direction of the dominant orientation:

   $θ = π 2 + tan − 1 2 J 12 J 22 − J 11 .$

   (16)
2. Bigün and Granlund⁶ described the anisotropy measure C as

   $C = λ 2 − λ 1 λ 2 + λ 1 ,$

   The domain kernel scale factors are set as

   $γ 2 = ( 1 + C ) and γ 1 = 1 / γ 2 .$

   (17)

   In constant neighbourhoods, where λ₁ + λ₂ = 0, C is set to 0 since there is no unique orientation. The aspect ratio increases with increasing anisotropy.

The structure tensor does not contain additional information than that present in the gradient, but it has the advantage that the matrix can be smoothed without cancellation effects in areas where gradients have opposite signs, since (∇I)(∇I)^T = (−∇I)(−∇I)^T. The oriented Gaussian kernel adapts to the data according to γ and θ obtained from the structure tensor and helps to smooth along edges. The oriented Gaussian kernel when combined with a range kernel (Fig. 6), results in the directional bilateral filter (DBF) given by

which is endowed with both edge- and direction-preserving capabilities.

In a practical scenario, since one does not have access to the original image x, we propose to use SURE to determine the optimal parameters ρ_d and ρ_r of the DBF. The differential of the filter output with respect to the noisy image is obtained as

Since the domain kernel weights are precomputed, the derivative of (18) with respect to y_p is approximated as

Using (19) and (20), we evaluate the divergence term and calculate $SURE ( x ˆ )$ in (9). The optimal parameters of the DBF are computed by minimizing SURE over a chosen range of possible values of the parameters ρ_d and ρ_r.

To further improve the performance of the directional bilateral filter, the optimal decay parameters of the Gaussian functions ρ_d and ρ_r are calculated over smaller patches of the image rather than the entire image. The use of patches is motivated by the observation that a fixed set of parameters ρ_d and ρ_r fixed globally may not be optimal for smaller patches of the same image. The image is divided into patches of size M × M, and the optimal decay parameters are calculated for each of the patches as described in Section IV. Typically, the patches are chosen to be one-eighth the size of the image and with 50% overlap between patches.

Our experiments involve examining whether the SURE follows the MSE for the bilateral filter, examining whether the estimate of the noise variance is close to the actual value, and qualitative and quantitative comparisons between the directional bilateral filter and the Gaussian bilateral filter. Before we show the results on microscopy images of biological specimens, we show the results on an image of silk fibers under different controlled noise conditions (Fig. 9(a)). The image contains important directionality information akin to that in F-actin fibers. An orientation map of the image is shown in Fig. 7. The orientation of the ellipses indicates the direction in which smoothing is performed. The aspect ratio indicates the local anisotropy – higher the anisotropy, higher the aspect ratio.

A noisy realization is obtained by adding zero-mean white Gaussian noise. The orientation and scaling parameters are computed from the noisy image using the structure tensor approach. The image is then denoised using the directional bilateral filter for different parameter settings ρ_d and ρ_r. In each case, the MSE and SURE were computed. The results are shown in Fig. 8. We observe that the SURE closely approximates the MSE. Further, the MSE curve is relatively flat around the optimal parameters for the DBF compared to the GBF. This makes it possible to use a coarser grid of parameters for the search.

The techniques are also compared with some of the state-of-the-art denoising methods such as SURE-LET,⁷ non-local means⁸ (NLM), and guided image filtering¹⁹ (guided). Image quality metrics such as PSNR and SSIM⁴³ are used for comparison. These approaches are termed full reference, meaning that a complete reference image is available for comparison. Since the technique is primarily developed for smoothing real microscopy images, we also use a no-reference or blind quality assessment metric called the BRISQUE²⁷ score, which is applicable in cases where the reference image is unavailable. The BRISQUE score typically has a value between 0 and 100 (0 represents the best quality, 100 the worst). For the guided image filtering technique, one iteration obtained by using the original image as the guide resulted in poor estimates. The performance was improved by using several iterations of the filter using the filtered output as the guide for the next iteration.²⁵ 

Quantitative comparisons are made based on the three quality metrics for various noise levels and are summarized in Table I. We observe that $σ ˆ$ is close to σ, that is, the estimate of the noise variance is close to the actual value. The PDBF has best performance over all noise levels in terms of PSNR. At σ = 10, it has the best performance in terms of all the metrics considered. However, as the noise level increases, NLM has better performance in terms of BRISQUE, and the guided image filter has better performance in terms of SSIM.

In the zoomed versions of the silk fibre image (Fig. 9), we observe that the edges appear sharper when filtered by the DBF (Fig. 9(e)) than the GBF (Fig. 9(d)). This is because of the directional smoothing. The optimal parameters for the GBF, DBF and PDBF are chosen by minimizing the SURE cost. For the PDBF, the optimal parameters are chosen locally, which further enhances its performance. This improvement is significant at higher noise levels. For all the real images acquired, the quantitative comparison is only in terms of BRISQUE as no reference images are available.

Cytoskeletal filamentous actin in the monolayer culture of the bovine pulmonary artery endothelial cells (BPAE line) was labeled with a fluorophore Alexa Fluor 488 (λ_exc ∣ λ_emi = 495 nm ∣ 519 nm) conjugated to phalloidin. Images were captured using a widefield as well as a confocal microscope.

The fluorophore was excited at 488 nm using an infinity-corrected 100x variable NA (1.3-0.6) objective lens (Olympus UplanFLN) for which the stoke-shifted emission was observed at 519 nm. The same objective lens was used to collect the fluorescence signal from the sample. Subsequently, the light was filtered by the band emission filter cube (Chroma, UBG triple-bandpass emission filter cube) to remove the scattered light. This was then focused onto a CCD camera (Jenoptik, ProGres MFcool).

Low NA (0.6) images were also obtained using a confocal microscope (Zeiss, LSM 700 AxioObserver) equipped with a low NA objective lens 40X (Zeiss, LD Plan Neofluor 0.6 NA) and a high NA oil immersion objective lens 63X (Zeiss, Plan Apochromat 1.4 NA) at different pinhole diameters (50 Airy units and 150 Airy units). A single objective lens was used to illuminate the sample as well to collect the fluorescence from the specimen. The fluorescence collected from the objective was efficiently filtered by the FixGate main beamsplitter. It was then directed through a fully automatic, high-precision pinhole that exclusively allows in-focal plane photons to fall on the detector. The light source is then moved to the next point on the specimen and the image is acquired again. The final image on the detector is the integrated effect of these points.

We assume the noise in these images to be additive white Gaussian (AWGN) and estimate the variance using (12). We filtered the image using the techniques shown in Table I. In Fig. 10 and Fig. 11, we show the denoised outputs for the top two BRISQUE scores as indicated in Table II. We considered several other widefield and confocal fluorescence microscopy images and the PDBF consistently resulted in lower BRISQUE scores than the other methods considered for comparison. Since the variances estimated in these images are low (σ < 10), we expect the proposed patch-based directional bilateral filter to have the best performance (cf. Table I).

Since the noise is assumed to be additive, the difference images (Fig. 10 and 11(d), and 11(e)) constitute the noise that has been suppressed. The noise removed by the SURE-LET and the NLM techniques (Fig. 10(e) and 11(e)) contains important directional information that has been removed while filtering. However, when the image is filtered by the PDBF, the difference images (cf. Fig. 10(d) and 11(d)) show no directionality information indicating that they have been preserved while smoothing. In Fig. 12, we show the advantage of filtering and denoising these images. The plot profiles along the yellow line (shown in the zoomed versions) show that the peaks that correspond to presence of filaments are well separated and can be resolved, which is a harder task to accomplish on the original noisy image.

We developed a directional bilateral filter for denoising and enhancing images obtained through fluorescence microscopy. This modified bilateral filter combines two edge-preservation techniques. The domain kernel incorporates orientation and anisotropy of image structures by means of a structure tensor and smoothes along the dominant orientations. When combined with the range kernel, the two kernels assist each other in edge preservation. We chose the optimal parameters of the directional bilateral filter by minimizing the SURE cost. The parameters that minimize SURE have been found to be nearly optimal in the MSE sense. We observe that choosing the parameters for smaller patches in an image improves the denoising performance. The proposed patch-based directional bilateral filter has better denoising performance than the Gaussian bilateral filter, SURE-LET, non-local means, and guided image filtering techniques in terms of PSNR. At low noise levels, it also had better performance in terms of SSIM and BRISQUE. We tested the filter on images containing cytoskeletal filamentous actin in the monolayer culture of the bovine pulmonary artery endothelial cells. Experiments show that directional bilateral filtering allows us to distinguish filaments easily. The directional bilateral filter efficiently suppresses noise while retaining important structural information making it a useful post-processing tool for microscopy images.

We would like to thank the reviewers and the editor for their suggestions and comments, which have helped improve the quality of the manuscript.