1.IntroductionThe 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.1 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 work2 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,3 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 production4 5 6 and its modulation by basic cellular factors such as, e.g., myosin isoforms,7 pH and regulatory proteins,8 or medically relevant substances such as, e.g., volatile anesthetics.9 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.10 A quantitative comparison of algorithms for tracking single fluorescent particles11 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. 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 filters12 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 Methods2.1.Anisotropic DiffusionThere 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,18 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 FilteringEvery 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. 2.1.2.Mathematical FormulationAnisotropic diffusion with a diffusion tensor evolves the initial image u under an evolution equation of type 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 :25 with convolution * , a Gaussian Gρ with standard deviation ρ, and uξ≔Gξ *u, a regularized version of u. The normalized eigenvectors ei 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.26 Using a diagonal matrix M with M ii=μi, Jρ can be written The diffusion tensor D uses the same eigenvectors ei. With the directional diffusivities λ1, λ2, λ3 and a diagonal matrix Λ with Λ ii=λi it becomes 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 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.21 So far the continuous equations for anisotropic diffusion have been described. We now proceed by describing their discretization.2.1.3.Discretization with Optimized FiltersEquation (1) can be solved using an Euler forward approximation for ∂u/∂t: where τ is the time step size and ui l denotes the approximation of u(x,t) in the voxel i at (diffusion) time l⋅τ. As defined above, Ai l denotes the spatial operator ∇⋅Di l∇. We use optimized separable first-order derivative filters to discretize Ai 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., h3 is the voxel volume) and 1-D smoothing kernels (e.g., [3,10,3]/16h) , denoted B in all other directions: 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(h2) 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.12We can further exploit the error term O(h2), which is always present. Numerical consistency does not define the behavior of Ai l for high spatiotemporal frequencies.12 Thus the properties of Ai 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 can be rewritten as: where χ is a scalar used as tuning parameter. If I1 and I2 were implemented identically the additional term χ(I2−I1) in Eq. (9) cancels out. As we want to modify the error behavior, I2 and I1 are implemented with different error behavior. Therefore, the first identity I1 is implemented by a scalar 1, thus without any error. The second identity I2 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 τ:
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. 2.2.In vitro Motility Assay2.2.1.Proteins and Experimental ChamberThe detailed description of protein isolation, solutions, and the preparation of the flow cell for the in vitro motility assay can be found in Ref. 2. In brief, rabbit skeletal muscle heavy meromyosin (100 μg/ml) was bound to the bottom of a microscopic flow chamber, consisting of a 22×50 mm 2 glass microscope slide and a 22×40 mm 2 precleaned coverslip coated with 0.1 nitrocellulose dissolved in amyl acetate. Rabbit skeletal muscle actin was prepared according to Pardee and Spudich28 with minor modifications as described in Ref. 2 and labeled with the fluorescent probe tetramethylrhodamin-phalloidin (R-415, Molecular Probes, OR, USA). Actin filaments (0.5 μg/ml) were added to the flow chamber and movement was initiated by adding a solution containing 2 mM ATP. All experiments were carried out at room temperature, ionic strength Γ/2=50 mM, and pH 7.4. 2.2.2.Fluorescence ImagingThe movement of the labeled actin filaments is observed with a high-resolution epifluorescence setup consisting of an inverted microscope (IX70, Olympus, Japan) equipped with epifluorescence illumination (Xe-light source) and a 100× objective (UPLANFL, 1.3 NA, oil, Olympus, Japan). The excitation bandpass filter is centered at 550 nm with a bandwidth of 10 nm, the beam splitter is at 560 nm, and the emission filter is centered at 580 nm with a bandwidth of 8 nm. The acquisition of the fluorescent image sequences is carried out with an intensified CCD camera (Luminescence Imager, Photonic Science, UK) with fiber optical coupling. The image sequences (512×512 pixel, 8 bit, 100 images) are digitized with 25 Hz temporal resolution using a PC-based frame grabber system (Meteor, Matrox, Canada). The analysis and further evaluation is carried out on a different PC-based computer system as described before. Additionally, the raw sequences are stored on video tape (AG-7355, Panasonic, Japan) for documentation. 3.ResultsAs 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 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.2 The addition of Gaussian noise to the test sequences generally leads to the broadening of the velocity distribution.2 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,29 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. 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 dimensions21 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 biweight27 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. 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. 4.DiscussionIn 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. AcknowledgmentsThis work was supported by grants from the Deutsche Forschungsgemeinschaft FOR240/3-1, Fi674/1-1 and the Bundesministerium fu¨r Bildung und Forschung BMBF 13N7871. REFERENCES
D. Uttenweiler
,
C. Veigel
,
R. Steubing
,
C. Goetz
,
S. Mann
,
H. Haussecker
,
B. Jaehne
, and
R. H. A. Fink
,
“Motion determination in actin filament fluorescence images with a spatio-temporal orientation analysis method,”
Biophys. J. , 78 2709
–2715
(2000). Google Scholar
S. J. Kron
and
J. A. Spudich
,
“Fluorescent actin filaments move on myosin fixed to a glass surface,”
Proc. Natl. Acad. Sci. U.S.A. , 83 6272
–6276
(1986). Google Scholar
M. Anson
,
M. A. Geeves
,
S. E. Kurzawa
, and
D. Manstein
,
“Myosin motors with artificial lever arms,”
EMBO J. , 15 6069
–6074
(1996). Google Scholar
G. Cuda
,
E. Pate
,
R. Cooke
, and
J. R. Sellers
,
“In vitro actin filament sliding velocities produced by mixtures of different types of myosin,”
Biophys. J. , 72 1767
–1769
(1997). Google Scholar
J. R. Sellers
,
G. Cuda
,
F. Wang
, and
E. Homsher
,
“Myosin-specific adaptations of the motility assay,”
Methods Cell Biol. , 39 23
–49
(1993). Google Scholar
M. Canepari
,
R. Rossi
,
M. A. Pellegrino
,
C. Reggiani
, and
R. Bottinelli
,
“Speeds of actin translocation in vitro by myosins extracted from single rat muscle fibres of different types,”
Exp. Physiol. , 84 803
–806
(1999). Google Scholar
A. M. Gordon
,
M. A. LaMadrid
,
Y. Chen
,
Z. Luo
, and
P. B. Chase
,
“Calcium regulation of skeletal muscle thin filament motility in vitro,”
Biophys. J. , 72 1295
–1307
(1997). Google Scholar
Y. Tsuda
,
T. Mashimo
,
I. Yoshiya
,
K. Kaseda
,
Y. Harada
, and
T. Yanagida
,
“Direct inhibition of the actomyosin motility by local anesthetics in vitro,”
Biophys. J. , 71 2733
–2741
(1996). Google Scholar
W. Hamelink
,
J. G. Zegers
,
B. W. Treijtel
, and
T. Blange
,
“Path reconstruction as a tool for actin filament speed determination in the in vitro motility assay,”
Anal. Biochem. , 273 12
–19
(1999). Google Scholar
M. K. Cheezum
,
W. F. Walker
, and
W. H. Guilford
,
“Quantitative comparison of algorithms for tracking single fluorescent particles,”
Biophys. J. , 81 2378
–2388
(2001). Google Scholar
H. Scharr
and
D. Uttenweiler
,
“3D anisotropic diffusion filtering for enhancing noisy actin filament fluorescence images,”
Lecture Notes in Computer Science , 2191 69
–75
(2001). Google Scholar
D. Uttenweiler
,
H. Scharr
,
C. Weber
,
B. Jaehne
, and
R. H. A. Fink
,
“Spatio-temporal anisotropic diffusion filtering of low S/N fluorescence image sequences,”
Biophys. J. , 80
(1), 504a
(2001). Google Scholar
V. Caselles
,
J. M. Morel
,
G. Sapiro
, and
A. Tannenbaum
,
“Special Issue on PDEs and Geom. Driven Diffusion,”
IEEE Trans. Image Process.,
(1998). Google Scholar
T. Preusser
and
M. Rumpf
,
“An adaptive finite element method for large scale image processing,”
Lecture Notes in Computer Science , 1682 223
–234
(1999). Google Scholar
G. H. Cottet
and
M. El Ayyadi
,
“A Volterra type model for image processing,”
IEEE Trans. Image Process. , 7 292
–303
(1998). Google Scholar
G. H. Cottet
and
L. Germain
,
“Image processing through reaction combined with nonlinear diffusion,”
Math. Comput. , 61 659
–673
(1993). Google Scholar
J. Weickert
and
H. Scharr
,
“A scheme for coherence-enhancing diffusion filtering with optimized rotation invariance,” Special Issue on PDE in Image Processing, Computer Vision and Computer Graphics,”
J. Visual Commun. Image Represent , 13 103
–118
(2002). Google Scholar
M. J. Black
,
G. Sapiro
,
D. Marimont
, and
D. Heeger
,
“Robust anisotropic diffusion,” Special Issue on PDEs and Geom. Driven Diffusion,”
IEEE Trans. Image Process. , 7 421
–432
(1998). Google Scholar
J. D. Pardee
and
J. A. Spudich
,
“Purification of muscle actin,”
Methods Cell Biol. , 24 271
–289
(1982). Google Scholar
|
CITATIONS
Cited by 26 scholarly publications and 2 patents.
Anisotropic diffusion
Anisotropic filtering
Diffusion
Image filtering
Luminescence
Signal to noise ratio
3D image processing