|Home | About | Journals | Submit | Contact Us | Français|
The interpretation and measurement of the architectural organization of mitochondria depend heavily upon the availability of good software tools for filtering, segmenting, extracting, measuring, and classifying the features of interest. Images of mitochondria contain many flow-like patterns and they are usually corrupted by large amounts of noise. Thus, it is necessary to enhance them by denoising and closing interrupted structures. We introduce a new approach based on anisotropic nonlinear diffusion and bilateral filtering for electron tomography of mitochondria. It allows noise removal and structure closure at certain scales, while preserving both the orientation and magnitude of discontinuities without the need for threshold switches. This technique facilitates image enhancement for subsequent segmentation, contour extraction, and improved visualization of the complex and intricate mitochondrial morphology. We perform the extraction of the structure-defining contours by employing a variational level set formulation. The propagating front for this approach is an approximate signed distance function which does not require expensive re-initialization. The behavior of the combined approach is tested for visualizing the structure of a HeLa cell mitochondrion and the results we obtain are very promising.
To date, it is firmly established that mitochondrial function plays an important role in the regulation of apoptosis (Green and Reed, 1998; Obeid et al., 2007). For instance, following a variety of cell death signals, mitochondria exhibit early alterations in function and morphologic changes, such as the opening of the permeability transition pore or mitochondrial megachannel (Frank et al., 2001; Zamzami et al., 2007). There is also strong evidence that defects in function may be related to many of the most common diseases of aging, such as Alzheimer dementia, Parkinson’s disease, type II diabetes mellitus, stroke, atherosclerotic heart disease, and cancer. This is founded on the observation that mitochondrial function undergoes measurable disturbance accompanied by drastic morphologic alterations in the presence of these multisystem diseases (Frey et al., 2006; Munnich and Rustin, 2001; Tandler et al., 2002).
Concurrent with the aforementioned conceptual advances there has been a significant increase in the types of tools available to study the correlation between mitochondrial structure and function. Along with the now classic methods for isolating mitochondria and assaying their biochemical properties, there are new and powerful methods for visualizing, monitoring, and perturbing mitochondrial function while assessing their genetic consequences (Marco et al., 2004; Pon and Schon, 2007). Electron tomography (ET) has allowed important progress in the understanding of mitochondrial structure. This imaging technique currently provides the highest three-dimensional (3D) resolution of the internal arrangement of mitochondria in thick sections (Perkins and Frey, 1999; Mannella et al., 1994). Nevertheless, the interpretation and measurement of the structural architecture of mitochondria depend heavily on the availability of good software tools for filtering, segmenting, extracting, measuring, and classifying the features of interest (Frey et al., 2002; Perkins et al., 1997).
This paper is organized as follows: section 2 presents an overview of anisotropic nonlinear diffusion models in image processing in general, and in electron microscopy in particular. The level set method is also presented briefly as it is applied to the extraction of contours in images. In section 3 we propose a new image smoothing and edge detection technique for electron tomography as an extension to the model proposed by Bazán and Blomgren (2007). This approach employs a combination of anisotropic nonlinear diffusion and bilateral filtering. In section 4 we exhibit the performance of the combined approach for visualizing the structure of a HeLa cell mitochondrion with very promising results. We end this paper with a summary and discussion in section 5.
In this section we present an overview of anisotropic nonlinear diffusion models in image processing and electron microscopy. The level set method is also presented briefly as it is applied to the extraction of contours in images. We only review here the works that serve as background to the model we propose in section 3. For an excellent and comprehensive survey of diffusion methods in image processing we refer the interested reader to the book by Weickert (Weickert, 1998) and the references therein. Two very good references for the level set method are the books by Osher and Fedkiw (2003) and by Sethian (1999).
Nonlinear diffusion is a very powerful image processing technique used for the reduction of noise and enhancement of structural features. It was first introduced to the image processing community by Perona and Malik (1990) as an attempt to overcome the shortcomings of linear diffusion processes, namely the blurring of edges and other localization problems. Their model accomplishes this by applying a process that reduces the diffusivity in areas of the image with higher likelihood of belonging to edges. This likelihood is measured by a function of the local gradient |u|. The model can be written as
for t > 0, on a closed domain Ω, with the observed image as initial condition u (x, 0) = u0 (x), and homogeneous Neumann boundary conditions g · u, n = 0, on the boundary Ω. Here, n denotes the outward normal to the domain’s boundary Ω, and g·u, n indicates the inner product ∫Ω (g·u) · n ds. In this model the diffusivity has to be such that g (|u|2) → 0 when |u| → ∞ and g (|u|2) → 1 when |u| → 0.
Notwithstanding the practical success of the Perona-Malik model, it presents some serious theoretical problems such as (i) ill-posedness (Nitzberg and Shiota, 1992; Weickert and Schnörr, 2000); (ii) non-uniqueness and instability (Catté et al., 1992; Kichenassamy, 1997); (iii) excessive dependence on numerical regularization (Benhamouda, 1994; Fröhlich and Weickert, 1994). The last observation motivated an enormous amount of research towards the incorporation of the regularization directly into the partial differential equation (PDE), to avoid too much implicit reliance on the numerical schemes. A variety of spatial, spatio-temporal, and temporal regularization procedures have been proposed over the years (Alvarez et al., 1992; Catté et al., 1992; Cottet and Germain, 1993; Weickert, 2001, 1994b, 1996b; Whitaker and Pizer, 1993). In subsection 2.2 we describe one of the variants to the Perona-Malik model that has been successfully used in electron microscopy, and in section 3 we propose a new model based on a combination of anisotropic nonlinear diffusion and bilateral filtering.
One way of introducing regularization to the Perona-Malik model is through anisotropic diffusion. Förstner and Gülch (1987) and Bigün and Granlund (1987) concurrently introduced the matrix field of the structure tensor for image processing, and it is the basis for today’s anisotropic diffusion models. The main idea behind these models is to construct the orthogonal system of eigenvectors v1, v2, of the diffusion tensor Dσ in such way that they will reveal the presence of edges, i.e., v1 || uσ (parallel) and v2 uσ (perpendicular). Then one chooses appropriate (corresponding) eigenvalues that will allow smoothing parallel to the edges and avoid doing so across them. The main advantage of anisotropic diffusion models over their isotropic counterparts is that they not only account for the modulus of the edge detector, but also its directional information.
where the vectors vi are the eigenvectors of the image’s structure tensor or its convolved version Jρ = Gρ *Jσ, where uσ = Gσ*u and Gσ, Gρ are Gaussian kernels of width σ, ρ, respectively. The parameters λi are functions of the eigenvalues, μ1 ≥ μ2 ≥ μ3, of the structure tensor Jσ (or Jρ). Together, the eigenvalues μi and the eigenvectors vi, characterize the local structural features of the image u, within a neighborhood of size O (ρ). Each eigenvalue μi reflects the variance of the gray level in the direction of the corresponding eigenvector vi, while each parameter λi controls the diffusion flux in the direction of vi and has to be chosen carefully.
Anisotropic nonlinear diffusion in electron microscopy was introduced by Frangakis and Hegerl (2001, 1999). Based on the works of Weickert (1998, 1999a,b), they chose the parameters λi to create a hybrid model that combines both edge enhancing diffusion (EED) and coherence enhancing diffusion (CED). EED is based on the directional information of the eigenvectors of the structure tensor Jσ, and its aim is to preserve and enhance edges. CED is based on the directional information of the eigenvectors of the convolved structure tensor Jρ, and is intended for improving flow-like structures and curvilinear continuities. To combine the advantages of EED and CED, the Frangakis-Hegerl model uses a switch based on comparing an ad hoc threshold parameter to the local relation between structure and noise (μ1 − μ3). The threshold parameter is based on the mean value of (μ1 − μ3) in a subvolume of the image containing only noise. EED is used when the difference (μ1 − μ3) is smaller than the threshold parameter. When it is larger, the model switches to CED for the enhancement of line-like structure patterns. In a separate publication, Frangakis et al. (2001) applied the hybrid model to two-dimensional (2D) and 3D electron tomography data and compared it with conventional methods as well as with wavelet transform filtering. They concluded that the model exhibits excellent performance at lower frequencies, achieving considerable improvement in the signal-to-noise-ratio (SNR) that greatly facilitated the posterior segmentation and visualization.
Fernández and Li (2003, 2005) proposed a variant to the model by Frangakis and Hegerl (2001, 1999) for ET filtering by anisotropic nonlinear diffusion, capable of reducing noise while preserving both planar and curvilinear structures. They provided their model with a background filtering mechanism that highlights the interesting biological structural features and a new criterion for stopping the iterative process. The Frangakis-Hegerl model diffuses unidirectionally along the direction of minimum change, v3, and efficiently enhances line-like structures (where μ1 ≈ μ2 μ3). It was argued by Fernández and Li (2003) that a significant number of structural features from biological specimens resemble plane-like structures at local scale. Therefore, they defined a set of metrics to discern whether the features are plane-like, line-like, or isotropic. The metrics they defined are:
which satisfy 0 ≤ Pi ≤ 1, ∀i and P1 + P2 + P3 = 1. In Eq. (3) μ1, μ2, and μ3 are the eigenvalues of the convolved structure tensor Jρ. These metrics are such that when P1 > P2 and P1 > P3, we have a plane-like structure; when P2 > P1 and P2 > P3, we have a line-like structure; and when P3 > P1 and P3 > P2, we have an isotropic structure. To achieve edge enhancing diffusion and linear and planar enhancing diffusion, the Fernández-Li model uses the parameters in Eqs. (4) and (5), respectively, for the diffusion tensor:
with user-defined free parameters α (regularization constant, typically set to 10−3) and C2 > 0 and C3 > 0. For the case of isotropic structure, the model employs what Fernández and Li (2005) call ‘background diffusion’ based on Gaussian smoothing.
The clever anisotropic nonlinear diffusion for the ET approach mentioned above also presents a practical drawback. Brox and Weickert (2002) argued that the linear structure tensor Jρ derived from Jσ by smoothing each component using a Gaussian kernel with standard deviation ρ, closes structures of a certain scale very well and removes the noise appropriately. However, it only preserves orientation discontinuities and does not preserve magnitude discontinuities, causing object boundaries to dislocate. Brox et al. (2005) argued that as soon as the orientation in the local neighborhood is not homogeneous, the local neighborhood induced by the Gaussian filter integrates ambiguous structure information. This information might not belong together and could lead to erroneous estimations. They proposed two alternatives to overcome this problem. The first solution involves the use of robust statistics for choosing one of the ambiguous orientations (van den Boomgaard and van de Weijer, 2002). The second solution is to adapt the neighborhood to the data by using the Kuwahara-Nagao operator (Bakker et al., 1999; Kuwahara et al., 1976; Nagao and Matsuyama, 1979). van den Boomgaard (2002) showed that the classic Kuwahara-Nagao operator can be regarded as a ‘macroscopic’ version of a PDE image evolution that combines linear diffusion with morphologic sharpening.
Based on the latter observation, Brox and Weickert (2002) proposed to address the aforementioned problem by replacing the Gaussian convolution by a discontinuity preserving diffusion method. This is obtained by considering the structure tensor Jσ as an initial matrix field that is evolved under the diffusion equation
∀ i, j, where the evolving matrix field uij (x, t) uses Jσ (x) as the initial condition for t = 0. The matrix D (A) = T (g (λi)) · TT is the diffusion tensor for A = T (λi) · TT. The latter represents a principal axis transformation of A with the eigenvalues λi as the elements of a diagonal matrix, diag (λi), and the normalized eigenvectors as the columns of the orthogonal matrix, T. For the diffusivity, g = 1 − exp (−c/(s/λ)8), for c > 0, and λ the contrast parameter. In Eq. (6), σ denotes the operator where Gaussian derivatives with standard deviation σ are employed. This approach tends to prevent boundary dislocations while keeping the desirable properties of the linear structure tensor. In section 3 we introduce a new approach based on anisotropic nonlinear diffusion and bilateral filtering for electron tomography. It allows noise removal and structure closure at certain scales, while preserving both the orientation and magnitude of discontinuities.
Osher and Sethian (1988) developed a framework relying on a PDE approach for modeling propagating interfaces. These methods have been applied to recover shapes of 2D and 3D objects from visual data, as shown by Malladi et al. (1996). This modeling scheme makes no a priori assumptions about the object’s shape and starts with an arbitrary function, propagating it in the direction normal to the curve along its gradient field with a certain speed, to recover shapes in the image. The level set formulation allows both forward and backward motion of the initial front through the creation of a higher dimensional function (x, t), where the initial position of the front is embedded as the zero level set. The evolution of the function (x, t) is then linked to the propagation of the front itself through a time-dependent initial value problem. Many implementations (Osher and Fedkiw, 2003) of the level set method utilize a zero (or initial) level set such that (x, t = 0) = ±d (x), where ±d (x) is the signed distance to the position of the front. Throughout the evolution of the front, in order to avoid the formation of shocks, very flat shapes, and/or very sharp shapes, a re-initialization process is often used periodically to restore a signed distance function. The process of re-initialization can be complicated and expensive. There is no simple way to determine how and when the level set function should be re-initialized to a signed distance function.
Li et al. (2005) presented a variational formulation whose propagating front is an approximate signed distance function yet does not require re-initialization. The variational energy functional consists of both an internal energy term that forces the level set function to be kept as an approximate signed distance function, and an external energy term that drives the zero level set toward the sought object contours in the image. The total energy functional is given by
The first term in the sum is the internal energy. It helps prohibit the deviation of from a signed distance function, where μ > 0 is the parameter controlling the effect of the penalizing the deviation. () is a metric that characterizes how close is to a signed distance function, e.g., t = sign (0) (1 − ||):
The second term in the sum of Eq. (7) is the external energy term that moves the zero level curve toward the object boundaries. Given an image u we can define the following edge indicator function where Gσ is the Gaussian kernel with standard deviation σ, g = (1 + |Gσ*u|2)−1
With this we can further specify our external energy term:
for constants: λ > 0, ν, and terms: g () = ∫Ω gδ () | |dx, and g () = ∫Ω gH (−) dx where δ (·) is the univariate Dirac function and H (·) is the Heaviside function. The energy term g() computes the length of the zero level curve of while g() is the weighted area on the interior of the zero level set and speeds up the curve evolution. The coefficient ν serves to control both the speed and direction of the curve propagation and should be chosen appropriately depending on the relative location of the initial contour to the object of interest. For an initial contour outside the object, ν should be a negative value so that the contours may shrink to the object boundary; whereas, a positive value should be chosen for ν if the initial contour is inside the object so that the contours might expand to the boundary.
The use of this energy functional completely eliminates the need for the expensive re-initialization as the evolution of the level set function is the gradient flow that minimizes the overall energy functional. The internal energy term maintains the level set function as an approximate signed distance function while the external energy term drives the propagation. The evolution equation is determined using calculus of variations to differentiate and setting its Gâteaux derivative equal to zero, yielding the steepest descent process for minimization of the functional :
In subsection 2.2 we discussed the application of anisotropic nonlinear diffusion in electron tomography. The approach used by Frangakis and Hegerl (2001, 1999) and Fernández and Li (2003, 2005) is based on a hybrid EED/CED denoising mechanism that performs very well on data containing low- to mid-frequency signal components. The technique greatly facilitates image enhancement for subsequent segmentation and improved visualization of complex biological specimens. In this section we propose a new image smoothing and edge detection technique for electron tomography as an extension to the model proposed by Bazán and Blomgren (2007). This approach employs a combination of anisotropic nonlinear diffusion and bilateral filtering. Jiang et al. (2003) introduced bilateral filtering for the removal of noise from biological electron microscopy data. They showed that bilateral filtering is a very effective mechanism for suppressing the noise in tomograms while preserving high resolution secondary structure features. To the best of our knowledge, Bajaj and Yu (2005); Bajaj et al. (2003) were the first ones to experiment with bilateral filtering coupled to an evolution driven anisotropic geometric diffusion PDE. They called their method ‘volumetric anisotropic diffusion model’ and it is based on a level set formulation that uses bilateral filtering as a pre-filtering step in order to obtain more precise curvature information. For the enhancement of 2D features, the model requires the selection of two free parameters that control the diffusion rate, and a threshold switch associated with the image. Our proposed model aims at incorporating the best of both approaches, anisotropic nonlinear diffusion and bilateral filtering, in a single computationally robust implementation that requires neither the choosing of parameters nor the setting of threshold switches. Furthermore, this approach permits preservation of both orientation discontinuities and magnitude discontinuities, and avoids structure dislocations. The model is equipped with the diffusion stopping criterion proposed by Bazán and Blomgren (2007), based on the second derivative of the correlation between the noisy image and the filtered image. (See Appendix A for details on this diffusion stopping criterion.)
The structure tensor of a 3D tomogram u is a symmetric positive semidefinite matrix, Jσ. This structure tensor is the most stable and reliable descriptor of the local structure of an image (Weickert, 1995). Similar to the approach in (Bazán and Blomgren, 2007), we propose using a refined estimate of the gradient of u at voxel x = (x, y, z) obtained by applying a bilateral filter in place of the Gaussian kernel. Bilateral filtering is a technique for smoothing images while preserving edges. The first application of this method is attributed to Aurich and Weule (1995) and it was subsequently rediscovered by Smith and Brady (1997) and Tomasi and Manduchi (1998). The basic idea underlying bilateral filtering is to combine domain and range filtering, thereby enforcing both geometric and photometric locality. The model can be expressed as
with the normalization constant
Typically, Gσd will be a spatial Gaussian that decreases the influence of distant pixels, while Gσr will be a range Gaussian that decreases the influence of pixels u (ξ) with intensity values that are very different from those of u (x), e.g.,
The parameters σd and σr dictate the amount of filtering applied in the domain and the range of the image, respectively.
The new structure tensor is therefore . The image’s local structural features can be determined by performing the eigen-analysis of the structure tensor Jbf where, as before, the eigenvalues provide the average contrast along the eigendirections, and the corresponding eigenvectors give the preferred local orientations. We will take advantage of the information provided by the structure tensor at each voxel x to devise a robust anisotropic structure enhancement model for 3D electron tomograms of mitochondria. Images of mitochondria contain many flow-like patterns and they are often perturbed by large amounts of noise (and may also include the artifacts related to the limited tilt range). Thus, it is necessary to denoise and enhance them by closing interrupted structures. To exploit the coherence and curvilinear continuity while connecting possible interrupted lines and surfaces, we average the structure tensor Jbf over a region by applying bilateral filtering in the form bf = Gbf * Jbf. The directional information is thereby averaged, although the structure of the region is still maintained by the discontinuity preserving bilateral filtering.
The diffusion tensor, Dbf 3×3, controls the smoothing across the 3D tomogram. We define the diffusion tensor as a function of the structure tensor bf,
where vi are the structure tensor’s eigenvectors. The eigenvalues of the diffusion tensor, λi, define the strength of the smoothing along the eigendirections, vi, and allow the application of different diffusion processes: (i) Linear diffusion or Gaussian smoothing is applied when λi = 1, ∀i; (ii) Nonlinear diffusion is applied if λi = g (|u|2 ), ∀i; (iii) Anisotropic nonlinear diffusion can be applied by setting the values λi so they reflect the image’s underlying local structure.
As mentioned in subsection 2.2, it is now common to use a hybrid approach that switches the diffusion process from EED to CED and vice versa, based on selected ad hoc thresholds. Switching to a third diffusion mode, Gaussian diffusion (GD), in areas where the image becomes predominantly isotropic (based on another ad hoc threshold) has also been suggested (Fernández and Li, 2003, 2005). We propose to use the anisotropic diffusion process where the model switches among the three modes, EED/CED/GD, automatically based on information extracted locally from the signal. The model can be regarded as ‘structure enhancing diffusion’ (SED), where the eigenvalues are defined as
The coherence measures (μ1 − μ2)2 and (μ1 − μ3)2 are computed based on the eigenvalues of the averaged structure tensor bf, and the parameters C2 and C3 act as thresholds such that structures where (μ1 − μ2)2 > C2 and (μ1 − μ3)2 > C3 are regarded as planar patterns, while structures where (μ1 − μ2)2 < C2 and (μ1 − μ3)2 > C3 are regarded as linear patterns. In practice, the logical ‘if μ1 = μ2 then’ and ‘if μ1 = μ3 then’ are unnecessary if we use exp (−Ci/((μ1 − μi)2 + ε))with i = 2, 3, for small ε.
In Eq. (15), g(|ubf|2)is a monotonically decreasing function such as Perona-Malik’s
with λ > 0 the typical contrast threshold parameter. There are several ways to set this parameter. Perona and Malik (1990) suggested using the idea presented by Canny (1986) and set λ as a percentile, p, of the image gradient magnitudes at each iteration (they recommended the value p = 90%.) A by-product of this approach is a decreasing λ, which has an stabilizing effect on the diffusion process (Mrázek, 2001).
The advantages of the proposed definitions for the diffusion tensor’s eigenvalues become evident when we perform a 2D analysis of the behavior of λ1 and λ2. Assume we are standing on an edge pixel of an image u [0, 1], as shown in Fig. 1. In this case, uy = 0 and 0 < |ux| ≤ 1, depending on the gray values in the two regions. The typical hybrid approach will use a chosen threshold to switch between EED and CED (see Eqs. (4) and (5)), and set the diffusion tensor’s eigenvalues to either λ1 = g (|u|2) or λ1 = α; and either λ2 = 1 or λ2 = α+(1 − α) exp (−C/μ1 − μ2)2). Assuming for example the values, α = 10−3 and C = 4.5 ×10−4, and considering that , we can plot λ1 and λ2 along uy = 0 and interpret the following (see Fig. 2):
Before applying the proposed SED model to a 3D electron tomogram of a mitochondrion, we apply the model to a more familiar 2D image and to a 2D slice of the electron tomogram of a HeLa cell mitochondrion below. The experiment demonstrates the excellent performance of the proposed model, where the structures of the images have been enhanced while most the noise has been removed. For the Clown image, we ran the four models BF, EED, CED and SED to the points of maximum similarity between the noise-free image and each of its filtered images. The similarity was measured by the correlation coefficient between the noise-free image and each of the filtered images, corr (f, u). In Fig. 3 we can observe that the SED model not only enhances the structures of the images, but also removes sufficient noise for the filtered image to move closer to the noise-free images than the other filtered images. For the 2D slice of the electron tomogram of a HeLa cell motochondrion, we ran the three models EED, CED, and SED to their stopping-points (see Appendix A). We can see in Fig. 4 that the SED model allows for good preservation of the edges and the smoothing of the background, without creating artificially enhanced details.
The electron tomogram employed in our experiments corresponds to a HeLa cell and was obtained from a 250 nm semi-thick section across a mitochondrion expressing cytochrome c-GFP. In the interest of research not discussed here, apoptosis was induced in the mitochondria with 100μM etoposide for 15 hours. The imaging occurred before the release of cytochrome c or loss of membrane potential allowing the maintenance of normal mitochondrion profiles; however, the treatment caused elongation of the crista junctions. The use of a semi-thick section is advantageous because it allows accurate depiction of the inner membrane topology of the mitochondrion.
The microscope used was the FEI Tecnai 12 Transmission Electron Microscope (TEM) with magnification set at 11000. The EM tomography single-tilt series 3D reconstruction was obtained from the semi-thick sample by progressively tilting the specimen and recording images using a Teitz 214 digital camera. The tilting was conducted in increments of 2 degrees over an angular range of ±60°; the angular range is limited by the geometry of the apparatus that holds the sample. Once the tilt-series was collected on the digital camera, the IMOD Software Suite (Kremer et al., 1996) was used to process the images and obtain the 3D reconstruction of the electron tomogram.
3D models are constructed using the electron tomogram. The 3D tomogram can be represented as a series of parallel sections, one pixel thick, of constant z. (Here, the image pixel size is 1.27 nm.) Constructing the models then requires the tracing of the membrane profiles of the outer membrane, inner membrane, and cristae structures in each of many parallel sections of the tomogram. These tracings form a stack of membrane contours that, when input to a computer display program, create a 3D model that can be rotated and viewed at any angle Frey et al. (2002). The tracing of the mitochondrial structures in the tomogram is currently done manually. However, by applying the proposed variational level set algorithm, the process is made less tracer-dependent.
After the 3D tomogram of the HeLa cell mitochondrion has been reconstructed, we apply the algorithm described in section 3 for the removal of noise and the enhancement of the structural features. This step is critical for the posterior segmentation and extraction of the structure-defining contours. The problem to solve is
where the diffusion tensor Dbf is given by Eq. (14), the eigenvectors vi for i = 1, 2, 3, are the eigenvectors of the structure tensor bf, and the eigenvalues λi for i = 1, 2, 3, are defined in Eq. (15). We apply the standard explicit finite difference scheme using central difference to approximate the spatial derivatives, and forward difference to approximate the time derivative. The condition for stability, assuming δx = 1, is given by δt = 1/6 (Weickert et al., 1998).
Fig. 5 and Fig. 6 show some results of the proposed approach applied to the 3D electron tomogram of the HeLa cell mitochondrion. The proposed approach achieves excellent noise reduction while preserving the salient edge features. In order to facilitate the extraction of the structures, we synthetically enhance the contrast by applying the confidence connected segmentation algorithm (Meier et al., 1997). In this context, this simple region-growing segmentation method produces sufficiently good results for the extraction stage, but more flexible methods such as the watershed technique (Volkmann, 2002) or the Chan and Vese (2001) algorithm can easily be substituted. In order to effectively remove unwanted small-scale features, including noise and extraneous cellular material, we label the components in the segmented Fig. 7(a), and keep the large-scale features corresponding to the mitochondrion of interest. A typical slice in the series of images processed for this paper contained over 300 distinct components. In all but one slice were the components of interest contained in the 3 largest labeled structures; in the ‘outlier’ case a small piece of the inner membrane was disconnected in the particular slice and was the 7th largest component. The result of filtering out the small scale features is shown in Fig. 7(b) and Fig. 7(c); in the latter we have ‘clipped’ external features which have become connected to the outer membrane due to noise and limitations in the imaging process. After segmentation, the features are extracted using the level set approach described in subsection 2.3.
We adopt the formulation by Li et al. (2005) presented in subsection 2.3 to extract the contours in the mitochondrion images. One significant advantage of this formulation is the liberty allowed in selecting the initial level set function. Traditionally, using level set methods requires the initial level set to be a signed distance function 0 so that re-initialization can be applied. However, with this need eliminated, a much simpler initial function may be defined:
given arbitrary Ω0, a subset in the image domain Ω where Ω0 is the set of points on the boundary of Ω0, and η > 0. For our implementation η = 4 is selected; however, most any constant would work. For the purposes here we use a Dirac function δ(x) in Eq. (10) that is slightly smoothed. We define the regularized Dirac function δε(x) as follows
and utilize ε = 1.5 for our implementation.
In implementing the proposed level set method, we carefully selected both our timestep τ and coefficient μ to be safely within the range required for stability ( as explained in Li et al. (2005)), τ = 5 and μ = .04. The level set functions were initialized as the function 0 defined in subsection 2.3 with η = 4 using selected regions Ω0. Fig. 8 shows the successful extraction of both the crista structures and outer membrane of the 718 × 763 pixel mitochondrion image. Identification of the interior structures was conducted separately from the identification of the outer membrane due to the required opposite direction of contour evolution for each. In order that the initial contour expand to identify the inner structures ν was chosen to be −25 and the evolution required 53 iterations, whereas for the second initial contour, used to shrink to identify the outer membrane, ν was chosen to be 10 and the evolution required 38 iterations. These selections allowed accurate visualization of the boundaries of interest. Note that the algorithm was run on the image in Fig. 5 (lower-left) and the resulting contours have been displayed on the original electron tomogram image slice. Also note that even though the SED model and the segmentation algorithm were run on the 3D tomogram, the contour extraction is done on a slice-by-slice basis.
We have presented a multi-stage approach for extracting the mitochondrial structures from electron tomograms. In particular, we apply the strategy to a 3D tomogram of a HeLa cell mitochondrion.
In the initial reconstruction, or noise reduction phase, we propose a structure enhancing anisotropic nonlinear diffusion strategy: the local structure tensor Jbf is formed from the gradient information of a bilaterally smoothed version of the current image. In order to close gaps in structures caused by imaging limitations, the local structure tensor is further smoothed with a bilateral filter, forming a smoothed version of the structure tensor, bf. This helps to preserve both orientation discontinuities and magnitude discontinuities. The eigenvectors vi, of the smoothed structure tensor form the basis for the diffusion tensor Dbf, where the eigenvalues are prescribed so that there is a smooth interpolation, rather than a hard threshold switching of the diffusion characteristics between image areas of differing structure properties.
After the noise reduction phase, we synthetically enhance the contrast by applying the confidence-connected segmentation algorithm. The connected components of the segmented images are labeled, and we discard the small scale components which correspond to either noise or unwanted extra-mitochondrial structures. Following which, structures are extracted using a variational level set formulation which includes a term that drives the level set function toward a signed distance function. This both simplifies the initialization of the algorithm and removes the need for re-initialization.
The extracted contours are visualized in Fig. 8. The left panel shows the successful result of extracting both the outer membrane and the inner structures; in the right panel we show a 3D-rendering of a ‘stack’ consisting of 25 extracted contours from 2D slices. The results are very encouraging. This computational approach is potentially much faster, and is more robust than the hand-tracing of structures.
The authors would like to thank Dr. Terry Frey and the San Diego State University Mitochondria Research Group for their input and for providing the images used in our computations. From this group, our special thanks go to Mei Sun and Mariam Ghochani. This work has been supported in part by NIH Roadmap Initiative award R90 DK07015.
Bazán and Blomgren (2007) proposed a new (very simple) diffusion-stopping criterion inspired by observation of the behavior of the correlation between the noise-free image and the filtered image, corr (f, u), and the correlation between the noisy image and the filtered image, corr (u0, u). Although the former measure is only available in experimental settings it helps validate the usefulness of the latter. The nonlinear diffusion process starts from the observed (noisy) image, u0 (x), and creates a set of filtered images, u (x, t), by gradually removing noise and details from scale to scale until, as t → ∞, the image converges to a constant value. During this process the cor relation between the noise-free image and the filtered image increases as the filtered image moves closer to the noise-free image. This behavior continues until it reaches a peak from where the measure decreases as the filtered image moves slowly towards a constant value. During the same process the correlation between the noisy image and the filtered image decreases gradually from a value of 1.0 (perfect correlation), to a constant value, as the filtered image becomes smoother (see Fig. A-1). By comparing both measures, we observe that as corr (f, u) reaches its maximum (the best possible reconstructed image), the curvature of corr (u0, u) changes sign. This suggests that a good stopping point of the diffusion process is where the second derivative of corr (u0, u) reaches a maximum. The performance of the stopping criterion can be observed in Fig. A-1 along with the reconstructed images of Lena and the Clown (Fig. A-2). We observe that the stopping criterion is almost optimal, allowing the diffusion process to stop near the point where the three filtering methods reach their best possible image reconstructions.
The proposed SED model, Eq. (17), is given by
The algorithm to implement the proposed SED model is as follows:
|δt 1/6||% time-step|
|σd 1||% domain Gaussian kernel’s width|
|σr 10−2||% range Gaussian kernel’s width|
|d 2||% orientation Gaussian kernel’s width|
|r 10−2||% magnitude Gaussian kernel’s width|
|u (x, 0) u0 (x), on Ω||% set initial condition|
|1 corr (u, u0)||% correlation measure|
|ubf Gbf * u||% convolve image with bilateral filter kernel|
|u [ux uy uz]T||% estimate gradients|
|ubf [(ubf )x (ubf )y (ubf )z]T||% estimate gradients|
|% magnitude of the gradients|
|% structure tensor|
|bf bf * Jbf||% convolve structure tensor with bilateral filter|
|vj, μj eigen (bf), j = 1, 2, 3||% eigenvectors and eigenvalues|
|λ prctile(|ubf |)||% contrast parameter|
|% diffusivity function|
|C2 sqtm(μ1, μ2)||% coherence parameter|
|C3 sqtm(μ1, μ3)||% coherence parameter|
|λ1 g||% prescribe eigenvalue|
|% prescribe eigenvalue|
|% prescribe eigenvalue|
|Dbf [v1 v2 v3] diag (λ1, λ2, λ3) [v1 v2 v3]T||% form diffusion tensor|
|Dbf · u, n 0, on Ω||% set boundary conditions|
|= · (Dbf · u)||% diffusion term|
|u u + δt||% evolve the image|
|i+1 corr (u, u0)||% update correlation measure|
|ĉi+1 2i+1||% stopping criterion|
|until ĉi+1 ≤ 0|
We are working on the implementation of a user-friendly package that will be distributed freely under an open source license. The current code is implemented in MatLab® and Parallel MatLab® and will be ported to a faster implementation using GPGPU. We hope to complete this phase of the project by the Fall of 2009.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.