1.

Introduction

The development of sophisticated fluorescence imaging techniques has led to an enormous growth in our understanding of cellular, subcellular, and molecular processes and function. Especially the high temporal and spatial resolution of modern fluorescence microscopic methods has led to their widespread use in biomedical, biotechnical, and biophysical research as a standard tool.¹ Common applications include regulation of cellular functions by second messengers, structural studies of cellular components, metabolic pathway analysis, and single molecule studies for the basic understanding of the complex molecular interactions involved in cellular processes.

Despite the advances in microscopic techniques, the analysis of fluorescence images is still very challenging due to several problems inherent to microscopic fluorescence images. Generally, when recording high spatially and temporally resolved dynamic processes as microscopic image sequences, the signal-to-noise ratio (SNR) is very low due to the limited amount of fluorescence photons available for detection. Therefore, image sequences of dynamic processes pose very high demands on automated methods for image analysis in general. Especially when quantifying motion, concepts have to be derived that yield the highest accuracy possible for this dynamic low light level data.

In previous work² we successfully demonstrated the use of an optical flow based orientation analysis method (3-D structure tensor method) for analyzing the velocity distribution of actin filament movement in the in vitro motility assay. In this experiment, originally devised by Kron and Spudich,³ fluorescently labeled actin filaments move over a surface of the immobilized motor protein myosin. It is routinely used for studying the interaction of the proteins actin and myosin during force production^{4 5 6} and its modulation by basic cellular factors such as, e.g., myosin isoforms,⁷ pH and regulatory proteins,⁸ or medically relevant substances such as, e.g., volatile anesthetics.⁹ This experiment is an ideal example for high spatially and temporally resolved fluorescence microscopic image sequence data, where motion has to be measured quantitatively and where the high amount of noise requires robust and precise algorithms for motion analysis. Figure 1 shows a typical example of in vitro motility assay image sequence data, where actin filament motion is visible as the displacement of rod-like objects. The zoomed image parts in the middle row visualize that the high amount of noise leads to a degradation of filaments (marked by arrows), which poses severe problems for standard threshold operations, as can be seen in the lower row of Figure 1. These large segmentation errors are one of the main reasons that many of the classical particle tracking approaches for motion analysis fail to produce reliable results in this kind of noisy image sequence data, as pointed out by several groups.¹⁰ A quantitative comparison of algorithms for tracking single fluorescent particles¹¹ has shown that different tracking approaches should be considered under the various experimental conditions, especially for different levels of noise. As mentioned before, the 3-D structure tensor method is ideally suited for the accurate and robust determination of the velocity distributions of moving actin filaments in this experiment, avoiding the various problems using particle tracking approaches. However, if particle properties are of interest, this pixel-based approach cannot be directly applied since the structure tensor method does not allow the direct analysis of object features.

Figure 1

Noise in microscopic fluorescence image sequences of single actin filaments. The top row shows three images of a fluorescence time series of actin filament motion in the in vitro motility assay. The motion of the fluorescently labeled actin filaments is visible as the displacement of rod-like objects. A magnification of the image area indicated by the rectangle is given in the middle row. The high amount of noise inherent to dynamic molecular fluorescence data leads to the degradation of filaments (marked by arrows). If a standard threshold algorithm is applied (lower row) this results in nonclosed objects (marked by arrows). The scale bar is 10 μm.

Therefore, we have now developed a method for increasing the signal-to-noise ratio and for an accurate object restoration in image sequences, which is based on a spatiotemporal anisotropic diffusion filtering method. It also uses the structure tensor for a 3-D analysis of the spatiotemporal image structure, which is subsequently used for a structure adopted image smoothing. In contrast to existing 2-D anisotropic diffusion filtering methods, the additional use of the temporal information of image sequences leads to a significant improvement in image denoising. Additionally, the use of accurate filters¹² and optimized discretization schemes, as described in Sec. 2, is crucial for an optimal performance of the image enhancement algorithm. As we will show, our method even allows the restoration of objects, which are degraded by the high amount of noise present in these molecular fluorescence images without morphological errors, a prerequisite for an accurate analysis of dynamic particle properties. In order to quantify the precision of the method under known conditions we have additionally applied it to computer-generated test data. Preliminary results of this work have been presented in conference talks.^{13 14}

2.

Materials and Methods

2.1.

Anisotropic Diffusion

There is a large amount of partial differential equation (PDE) methods for image processing described in literature,^{15 16 17} including anisotropic diffusion. Numerical techniques include finite elements,¹⁸ lattice Boltzmann techniques (Jawerth et al.; in Ref. 15), and finite difference schemes.^{19 20 21 22} The numerical scheme used here was previously presented for the 2-D case.^{23 24} Its most prominent property is excellent directional behavior due to the use of optimal derivation filters. As motion in image sequences is equivalent to orientation in a spatiotemporal image stack, directional accuracy is crucial for denoising without changing the speed of moving object. Thus we have extended this scheme to three dimensions and further enhanced it by additionally considering the error term (see below).

2.1.1.

Principle of Anisotropic Diffusion Filtering

Every dimension of an image sequence will be treated as a spatial dimension, thus a 2-D image sequence is treated as a 3-D data set. The original image will be changed by applying anisotropic diffusion, which is expressed by a diffusion time t.

It is well known in image processing that isotropic smoothing (or blurring) with a Gaussian or binomial filter corresponds to linear diffusion ∂u/∂t=DΔu, where u(x,t) denotes the original image, t is the diffusion time, Δ is the Laplace operator, and D corresponds to the diffusion coefficient. Thus for the scope of this paper anisotropic nonlinear diffusion can be considered as anisotropic smoothing by a deformed Gaussian as illustrated in Figure 2. Therefore smoothing and diffusion will be used synonymously in the remainder of the paper. It should be noted that when using anisotropic diffusion one also needs to carefully consider further properties of this method as, e.g., influences of the so-called shock term or backward diffusion. A basic assumption for the successful use of the method is that “signal” consists of oriented structures and “noise” is uncorrelated and thus has no preferred orientation. Consequently, fast damping of uncorrelated structures while preserving oriented structures will result in efficient denoising. We achieve this aim by large diffusion along oriented structures and almost no diffusion perpendicular to them.

Figure 2

Principle of anisotropic diffusion filtering. Image smoothing is considered as a nonlinear anisotropic diffusion process with structure adopted diffusion coefficients. This description is equivalent to image smoothing with Gaussian filters, where the smoothing kernel is structure adopted. The adaptation is such that smoothing is isotropic for areas without object structures and for edges and corners smoothing is primarily parallel to the direction of the structure and minimal perpendicular to it. Note: for simplicity the principle is shown for the two spatial image dimensions, whereas the method presented in this paper is a full 3-D anisotropic diffusion filter, with the two spatial image coordinates and the time t as the third image coordinate.

2.1.2.

Mathematical Formulation

Anisotropic diffusion with a diffusion tensor evolves the initial image u under an evolution equation of type

Eq. (1)

$$\frac{\partial u}{\partial t}=\nabla \cdot (\mathit{D}\nabla \mathit{u}),$$

with the evolving image u(x,t), diffusion time t, 3-D derivation vector ∇=(∂_x1,∂_x2,∂_x3), and diffusion tensor D, a positive definite, symmetric 3×3 matrix. It is adapted to the local image structure measured by the structure tensor J :²⁵

Eq. (2)

$${\mathit{J}}_{\rho }(\nabla {\mathit{u}}_{\xi })={{\mathit{G}}_{\rho }}^{*}(\nabla {\mathit{u}}_{\xi }\nabla {{\mathit{u}}_{\xi }}^{\mathit{T}}),$$

with convolution ^* , a Gaussian G_ρ with standard deviation ρ, and u_ξ≔G_ξ ^*u, a regularized version of u. The normalized eigenvectors e_i of J_ρ give the preferred local orientations, and its eigenvalues μ_i give the local contrast along these directions. The structure tensor is highly robust under isotropic additive Gaussian noise.²⁶ Using a diagonal matrix M with M _ii=μ_i, J_ρ can be written

Eq. (3)

$${\mathit{J}}_{\rho }=({\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3)\mathbf{M}{({\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3)}^T.$$

The diffusion tensor D uses the same eigenvectors e_i. With the directional diffusivities λ₁, λ₂, λ₃ and a diagonal matrix Λ with Λ _ii=λ_i it becomes

Eq. (4)

$$D=({\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3)\mathbf{\Lambda }{({\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3)}^T.$$

The directional diffusivities λ_i determine the behavior of the diffusion. They shall be high for low values of μ_i and vice versa. In our application we use

Eq. (5)

$${\lambda }_{\mathit{i}}≔\{\begin{array}{cc}1\hfill & \text{if}{\mu }_{\mathit{i}}\le \sigma \hfill \\ 1-(1-\mathit{c})\text{exp}\left[\frac{-\mathit{d}}{{({\mu }_{\mathit{i}}-\sigma )}^2}\right]\hfill & \text{else}\hfill \end{array},$$

where c∈]0,1], d>0, and σ>0 corresponds to the global absolute noise level. The condition number of D is bounded by 1/c. Instead of this choice other functions can be used. In comparative tests (see Sec. 3) we also use nonlinear isotropic diffusion with Tuckey’s biweight in two and three dimensional (see, e.g., Ref. 27) and edge-enhancing diffusion in two dimensions.²¹ So far the continuous equations for anisotropic diffusion have been described. We now proceed by describing their discretization.

2.1.3.

Discretization with Optimized Filters

Equation (1) can be solved using an Euler forward approximation for ∂u/∂t:

$$\frac{u_i^{l+1}-u_i^l}{\tau }=\nabla \cdot (D\nabla u_i^l)$$

Eq. (6)

$$\begin{array}{ccc}\hfill \iff {\mathit{u}}_{\mathit{i}}^{\mathit{l}+1}& \hfill =\hfill & [1+\tau ({\partial }_{x1},{\partial }_{x2},{\partial }_{x3})D{({\partial }_{x1},{\partial }_{x2},{\partial }_{x3})}^T]{\mathit{u}}_{\mathit{i}}^{\mathit{l}}\hfill \\ \hfill & \hfill =:\hfill & (1+\tau {\mathit{A}}_{\mathit{i}}^{\mathit{l}}){\mathit{u}}_{\mathit{i}}^{\mathit{l}},\hfill \end{array}$$

where τ is the time step size and u_i ^l denotes the approximation of u(x,t) in the voxel i at (diffusion) time l⋅τ. As defined above, A_i ^l denotes the spatial operator ∇⋅D_i ^l∇. We use optimized separable first-order derivative filters to discretize A_i ^l (see Ref. 12 for the 2-D case). They are composed of a common 1-D derivative stencil (e.g., [−1,0,1]/2h) , denoted D, where h is the spatial discretization step (i.e., h³ is the voxel volume) and 1-D smoothing kernels (e.g., [3,10,3]/16h) , denoted B in all other directions:

Eq. (7)

$${\partial }_{{\mathit{x}}_{\mathit{i}}}={\mathrm{D}}_{{\mathit{x}}_{\mathit{i}}}*{\mathrm{B}}_{{\mathit{x}}_{\mathit{j}}}*{\mathrm{B}}_{{\mathit{x}}_{\mathit{k}}}+O(h^2),$$

where {i,j,k} is a permutation of {1,2,3}, lower indices at the brackets give the direction of the kernel,  *  is a convolution, and O(h²) is the numerical discretization error term of order 2 vanishing for h→0. These filters approximate rotation invariance significantly better than related popular stencils by an order of magnitude.¹²

We can further exploit the error term O(h²), which is always present. Numerical consistency does not define the behavior of A_i ^l for high spatiotemporal frequencies.¹² Thus the properties of A_i ^l can be influenced there without altering numerical consistency order in the following way: we add and subtract the identity operator I in Eq. (6), and therefore

Eq. (8)

$${\mathit{u}}_{\mathit{i}}^{\mathit{l}+1}=(1+\tau {\mathit{A}}_{\mathit{i}}^{\mathit{l}}){\mathit{u}}_{\mathit{i}}^{\mathit{l}}$$

can be rewritten as:

Eq. (9)

$${\mathit{u}}_{\mathit{i}}^{\mathit{l}+1}=\{1+\tau [{\mathit{A}}_{\mathit{i}}^{\mathit{l}}+\chi ({\mathit{I}}_2-{\mathit{I}}_1)]\}{\mathit{u}}_{\mathit{i}}^{\mathit{l}},$$

where χ is a scalar used as tuning parameter. If I₁ and I₂ were implemented identically the additional term χ(I₂−I₁) in Eq. (9) cancels out. As we want to modify the error behavior, I₂ and I₁ are implemented with different error behavior. Therefore, the first identity I₁ is implemented by a scalar 1, thus without any error. The second identity I₂ is implemented with a discretization error of order 4 as separable convolution kernel [−1,4,10,4,−1]/16h applied to the three spatiotemporal directions x,y,t. * Consequently we solely get a modification in the error term and the update rule in Eq. (9).

The following algorithm is repeated for every time step τ:

1. (1) Calculate the structure tensor J [Eq. (2)].

2. (2) Get the diffusion tensor D by J [Eqs. (3) and (4)].

3. (3) Calculate the flux j_i≔∑_m=1 ⁿD_i,m∂_x,mu, ∀i∈{1,…,n}.

4. (4) Calculate A_i ^lu_i ^l=∇⋅(D∇u_i ^l) by A_i ^lu_i ^l=∑_m=1 ⁿ∂_xm j_m.

5. (5) Calculate I₂u_i ^l by convolving u with [−1,4,10,4,−1]/16h in x,y, and t direction.

6. (6) Update in an explicit way by u_i ^l+1=[1+τA_i ^l+χ(I₂−1)]u_i ^l.

The iteration number is N and the total diffusion time T is T=τN.

The anisotropic diffusion algorithm is implemented using the image processing software Heurisko 4.08 (AEON, Hanau, Germany) on a PC-based system (Intel Pentium III, 700 MHz, 1 GB RAM). Computation time for the anisotropic diffusion image denoising of an image sequence (512×512, 75 frames) is 208 s per iteration. Usually between 5 and 50 iterations suffice for most denoising tasks.

3.

Results

As shown in Figure 3 we have first applied the method to computer-generated test sequences in order to characterize the performance of the anisotropic diffusion filter method under known conditions. The test sequences are similar to the ones used in Ref. 2. In brief, rod-like objects with a gray value of 200 move against a black background with gray level 0. The objects move with defined displacements of 1 pixel/frame for the objects moving in x and y directions and in approximate circles and with $\sqrt{2}\text{ pixel/frame}$ for the objects moving in a 45-deg angle./ These two velocity populations are visible in the velocity histogram for the test data without noise in the left column of Figure 3. All histograms are obtained with the structure tensor method.² The addition of Gaussian noise to the test sequences generally leads to the broadening of the velocity distribution.² The level of noise added to the test data in the middle column of Figure 3 has been chosen such that it is larger than the limit where the structure tensor and classical particle tracking algorithms yield reliable results for the velocity distribution. The addition of Gaussian noise with a standard deviation of 70 gray values to the test sequences (middle column in Figure 3) leads to a velocity distribution, where both velocities populations cannot be separated by automated algorithms. Applying the spatiotemporal anisotropic diffusion filtering to the noisy image sequence yields the denoised image data shown in the right panels. The signal-to-noise level is significantly improved and objects can be easily detected without morphological changes in object shape induced by the algorithm. Defining the signal-to-noise ratio as SNR=20^* log[(g−b)/Δg], where g is the mean gray value of an object, b the mean background noise level, and Δg the standard deviation of the gray values of an object,²⁹ the signal-to-noise-ratio for the filament-like structures in the original noisy image sequence is computed as SNR=2.0 in comparison to SNR=7.6 for the data obtained with the anisotropic diffusion filtering. Hence, an improvement in signal-to-noise ratio by a factor 3.8 is achieved.

Figure 3

Test pattern used to analyze the accuracy of the anisotropic diffusion filtering method. The rod-like objects (gray value 200) move in the x, y direction, in a 45-deg angle, and in approximated circles, respectively, versus a black background (gray value zero). The test pattern with no noise is shown on the left and the test pattern with added Gaussian noise with a standard deviation of 70 gray values is shown in the middle panel. Application of the 3-D anisotropic diffusion filtering to the middle panel resulted in the panel shown on the right. In the middle row the magnifications of an image area, which is indicated by the rectangle in the original image, clearly show that the anisotropic diffusion filtering successfully restores the objects degraded by the added noise, while simultaneously eliminating the background noise. The respective velocity histograms as obtained with the structure tensor method² for the moving test objects are given below each panel. The restoration of the two velocity populations is successfully achieved by the 3-D anisotropic diffusion filtering and furthermore the peak velocities are very precisely restored.

The noise reduction is also reflected in the velocity histogram of the processed data, where the artifacts in velocity determination originating from noise (velocities<1 pixel/frame) are drastically reduced. Moreover, the anisotropic diffusion filtering successfully restores the two velocity populations present in the test data. Most importantly, the mean velocities of both filament movements, as derived from the mean μ of a Gaussian data fit applied to both peaks in the histogram, are very precisely restored (μ=1.04 pixel/frame with a broadness of the distribution of 0.12 pixel/frame, as obtained from twice the standard deviation σ, and μ=1.42 pixel/frame with a broadness of 0.22 pixel/frame). The error in peak velocity is below 5 , which is remarkable, since the two filament populations could not be separated in the original noisy data. This result proves that the spatiotemporal anisotropic diffusion filtering does not introduce velocity artifacts into the image sequences, a prerequisite for the use of a filter method in quantitative analysis of motion in image sequences.

The power of the 3-D anisotropic diffusion filter method in denoising fluorescence microscopic image sequences of single actin filaments is shown in Figure 4. Additionally a comparison with commonly used binomial and diffusion filtering schemes is also given (all schemes were implemented with optimized filters and accurate discretization). Row A shows the unprocessed original fluorescence image data demonstrating the high amount of noise resulting in filament degradation. Applying a standard 5×5 binomial filter results in the images shown in row B. The SNR is still very low. Using an edge-enhancing diffusion scheme in two dimensions²¹ leads to closed object structure as can be seen in row C. However, the filament morphology is significantly altered by this method, which introduces arbitrary errors when particle properties are of interest. The result of nonlinear isotropic diffusion with Tuckey's biweight²⁷ in two dimensions is shown in row D and the same scheme implemented in three dimensions yields the images in row E. The 2-D scheme is not able to restore closed object structures. Therefore, it does not significantly improve the ability to accurately perform an object segmentation in these image sequences. The 3-D scheme with the additional temporal information used for smoothing significantly improves the results and restores closed object structures except for heavy filament degradation. Finally, the anisotropic diffusion filter scheme in three dimensions as shown in row F yields closed objects without introducing morphological changes in filament shape for almost all filament degradations in these image sequences.

Figure 4

Comparison of various filtering methods. (A) Unprocessed original fluorescence microscopic data of actin filament movement in the in vitro motility assay. (B) Smoothing with a 5×5 binomial filter. (C) Edge-enhancing diffusion filtering in two dimensions. (D) Isotropic nonlinear smoothing with a spatiotemporal diffusion filter using Tuckey‘s biweight in two dimensions and (E) in three dimensions. (F) Spatiotemporal anisotropic diffusion filtering with the new scheme yields the best results for reducing the noise and restoring the morphology of objects.

In Figure 5 a more detailed analysis of the denoising schemes is exemplarily given for the filament marked by a solid arrow in Figure 4. The gray level intensities are plotted perpendicular to the filament axis (panels A to F) and along the filament axis (panels A2 to F2). It can be clearly seen that only the anisotropic diffusion scheme (panel F2) restores a homogenous gray value profile along the filament axis, whereas all other methods retain some filament discontinuities, making it difficult to set automatic threshold levels for an accurate and reliable object segmentation. As previously described for the computer-generated test data in Figure 3, the calculation of the signal-to-noise ratio for the original noisy motility assay image data and the data processed with the anisotropic diffusion filter method (panel F2) again yields a significant improvement of the signal-to-noise ratio by a factor 3.8. For the edge-enhancing diffusion algorithm in panels C and C2 to restore a homogenous gray value profile it is necessary to set the smoothing parameters such that it introduces heavy filament morphological errors due to an over-smoothing. The overall tendency of the various methods to blur structures can be quantified using the perpendicular intensity cuts. The values of the full widths at half maximum (FWHM) obtained from Gaussian fits applied to the data are as follows. A: FWHM 0.34 μm (this value corresponds to the FWHM of the microscopic point spread function), B: FWHM 0.50 μm, C: FWHM 0.92 μm, D: FWHM 0.93 μm, E: FWHM 0.45 μm, and F: FWHM 0.57 μm. This data is a representative example for the fact that the new anisotropic diffusion filter scheme performs best in restoring a homogenous gray value profile along the filament axis with a moderate tendency to blur edges and corners. It is also best suited to detect structures close to the noise level as, e.g., the filament marked by the dashed arrow in Figure 4. Thus, the data obtained with the 3-D anisotropic diffusion filter scheme allows a reliable segmentation of objects and the analysis of filament properties without significant artifacts due to noise or morphological changes introduced by the denoising scheme.

Figure 5

Exemplary comparison of filament morphology obtained with the various filtering methods. The filament morphology is quantified by plotting the gray values perpendicular and along the filament axis of the filament marked by a solid arrow in Figure 4. The perpendicular intensity cuts are given in panels A to F, corresponding to the nomenclature in Figure 4. The longitudinal intensity cuts are depicted by panels A2 to F2, respectively. For a better comparison the mean background noise level has been added to the plots (solid line) as well as the mean gray value level 40 above background noise (dashed line).

4.

Discussion

In summary we have presented a new method for the enhancement of noisy fluorescence microscopic image sequences. The method is based on a 3-D anisotropic diffusion filtering scheme, with the two image coordinates and the time t of the image sequence as the third image coordinate. Using fluorescence microscopic image sequences of moving actin filaments in the in vitro motility assay, we could demonstrate that this approach successfully increases the SNR, allowing a quantitative analysis of dynamic object properties, which is not possible with the raw data or when denoising the data with commonly used filtering schemes. The 3-D anisotropic diffusion scheme restores objects in noisy fluorescence image sequences without changing the morphology of the objects of interest. Furthermore, the precision of the scheme has been validated on computer-generated test sequences. The accuracy and the restoration power of the 3-D anisotropic diffusion scheme is mainly based on (1) the novel use of optimized adaptable filters and accurate discretization schemes, which is particularly crucial for filtering schemes in 3 dimensions, and (2) the additional use of the temporal information inherent in image sequences, which drastically improves the power of reducing noise without blurring artifacts, resulting in an accurate object restoration. Also the temporal characteristics are not altered by this method (e.g., velocities, event durations) since the smoothing kernel is also structure adopted in the temporal domain. It should be noted that an accurate implementation of the 3-D schemes is absolutely necessary for precise, stable, and robust performance. The high amount of computing power which is necessary for the 3-D anisotropic diffusion scheme should not be a limiting factor in the future, as computer power will continue to increase, leading to computing times in the minute range even on a single PC for typical image sequence sizes in the near future. In general, this method is independent of object shapes and therefore it is also applicable to other molecular and cellular studies such as, e.g., the analysis of the movement of vesicles and other spherical objects. Also, the extension to multidimensional data is generally possible with this anisotropic diffusion filtering approach.

In the context of the growing importance of fluorescence microscopic techniques in the entire field of life science, the method presented here significantly improves and extends these methods for dynamic low light level applications. It is of general use for all kinds of microscopic fluorescence image sequence data, where the high amount of noise especially in high spatially and temporally resolved image sequences has so far limited the quantitative use and accurate analysis of the data.