Feature Extraction of 3T3 Fibroblast Microtubule Based on Discrete Wavelet Transform and Lucy-Richardson Deconvolution Methods.

Accompanied by the increasing requirements of the probing micro/nanoscopic structures of biological samples, various image-processing algorithms have been developed for visualization or to facilitate data analysis. However, it remains challenging to enhance both the signal-to-noise ratio and image resolution using a single algorithm. In this investigation, we propose a composite image processing method by combining discrete wavelet transform (DWT) and the Lucy-Richardson (LR) deconvolution method, termed the DWDC method. Our results demonstrate that the signal-to-noise ratio and resolution of live cells' microtubule networks are considerably improved, allowing the recognition of features as small as 120 nm. The method shows robustness in processing the high-noise images of filament-like biological structures, e.g., the cytoskeleton networks captured by fluorescent microscopes.

1. Introduction

1.1. Research Background

In recent years, as key components, microdevices have been widely used in the development of biological and biomedical techniques for manipulating, mixing, separating, and sensing biological targets. At the same time, increasingly stringent requirements have been placed on imaging techniques involving high signal-to-noise ratio and structural resolution. Currently, the most commonly used devices for the high-resolution imaging of biological or biomedical targets include confocal microscopes [1], stimulated emission depletion (STED) microscopes [2], and structured light illumination microscopes (SIM) [3] etc. Furthermore, many algorithms have been developed to improve the spatial resolution and signal-to-noise ratio (SNR) of biological images, including degenerate-model-based algorithms (e.g., deconvolution [4,5,6,7,8]), mathematical transformation-based algorithms (e.g., spectrum analysis [9,10], DWT analysis [11,12,13,14,15,16]), and machine-learning-based algorithms (e.g., deep learning [17,18,19]). However, most of these algorithms focus on a single task, e.g., inhibiting noise, identifying structure contours, or improving resolution. Furthermore, these methods normally require the target images to be clear, with relatively low levels of noise.

Nevertheless, many fundamental and representative biological structures are small and irregular and suffer from significant noise backgrounds during imaging. For example, the microtubule of the fibroblast is a fundamental cell structure and plays an important role in cellular responses to external stimuli. It has filament-like structures with a width of ~25 nm [20] and forms a densely packed network in nature [21]. These factors make it difficult to distinguish a single microtubule filament and track its dynamics during various biological processes, especially against noisy backgrounds, in which fluorescence signals due to emitted background light and autofluorescence originate from the areas above and below the focal plane can decrease the SNR of image. For isotropic or quasi-isotropic features, e.g., round-shaped and nanometer-sized exosomes, a deconvolution-based algorithm can effectively improve the structural resolution [22,23]. For densely packed networks (e.g., the microtubule), researchers are still paying more attention on high-performance image-processing methods [24], to achieve higher resolution and SNR.

1.2. Previous Works

Deconvolution methods, including the Lucy–Richardson (LR) algorithm [25,26,27,28], the fast thresholded Landweber (FTL) algorithm, the generalized expectation maximization (GEM) algorithm, etc., are commonly used in improving the image resolution and quality. The LR method was widely used in image processing accompanied by STED microscopy. For instance, in 2006, Willig et al. studied green fluorescent protein (GFP)-labeled viruses with STED microscopy and the LR method, achieving a lateral resolution of about 70 nm [29]. In 2007, Hell used the LR method to increase the spatial resolution of STED microscopy to 20–30 nm [30]. The GEM algorithm was advanced in 2006 by Bioucas-Dias et al. to process macroscopic images [31]. Although the authors found that the GEM method could improve image quality, the algorithm only compares the SNR before and after processing, which cannot ensure the original image intensity distribution before and after image processing. The FTL algorithm is a fast variational deconvolution algorithm that minimizes a quadratic data term. Vonesch et al. used FTL to process confocal images of a neuron cell [32]. They found that the FTL algorithm could achieve an eight-decibel improvement in ten iterations, with an insignificant increase in the image SNR. However, deconvolution methods by themselves may lead to over-processing and spurious images, especially in images with poor SNR.

By contrast, the wavelet method has been commonly applied for denoising (e.g., the expectation maximization (EM) algorithm [33,34]) and the extraction of featured structures through scales [35]. For instance, the EM algorithm utilizes both wavelet transform and fast Fourier-transform to improve the SNR of images. It can increase the SNR of a macroscopic image from 3 dB to ~7 dB after 8 to 10 iterations.

2. Methodological Principles and Process

2.1. Target of Image Process

In this investigation, we developed a new algorithm named as DWDC method, which combines DWT and Lucy–Richardson deconvolution. With this method, the spatial resolution of a typical biological image (high noise, blurred, and unclear) can be increased with improved SNR, and the features of filament-like structures can be extracted.

Figure 1a shows a confocal fluorescence image of 3T3 fibroblast microtubule networks, which were taken using Nikon A1 microscope and Olympus 100X oil immersion lens (NA 1.4) (Nikon Corporation, Tokyo, Japan). The excitation light wavelength is 640 nm, and the emission peak is around 674 nm for the SiR-Tubulin dye (Cytoskeleton, Inc. Denver, CO, USA). Each fluorescence image has 512 × 512 pixels, with a dot pitch of 0.25 µm. The fluorescence image’s bits per pixel (BPP) is 16 and the contrast ratio (CR) is 6000:1. For 3D reconstruction, twenty images were captured by z-stacking, with a 1-micrometer vertical interval. It is obvious that the branch of microtubule structures is highly contaminated by noise and the structures are clearly bold (Figure 1b,c). We can hardly distinguish the topological structures due to the fragmental distribution of fluorescent intensity. Structural features reflecting cell–cell interactions are indistinguishable from the figures. Since the featured scale of biological samples is comparable to or even smaller than the image pixel pitch (250 nm), the actual distribution of biological structures can be affected by the noise distribution. It is necessary to remove as much noise as possible before analyzing biological activities.

Figure 1

Original confocal image of 3T3 cell microtubule. (a) Image has 512 × 512 pixels, with a dot pitch of 250 nm. (b) The zoomed-in image of the white box. (c) The three-dimensional (3D) reconstruction of (b). The white scale bar represents 10 µm (for further details on the 3D structure of Figure 1 (c), see Supplementary video S1).

2.2. Methods and Process

An optical image is a convolution of an object with the point spread function (PSF) of an optical system [36]. If is the matrix of the image, where is the PSF of the optical system, is the light distribution according to the object and is the measurement noise of the optical system. If the size of the PSF is larger than the size of the mesostructure of the actual object, the imaging result has an insufficient spatial resolution to reveal the detail of the original object. Accordingly, the image after the optical system is blurred relative to the actual object.

The DWDC method advanced in this investigation utilizes both LR and DWT, as presented in Figure 2. Firstly, the image was processed using Gaussian interpolation and threshold analysis. DWT was then applied to suppress noise level and extract characteristic microtubule structures on the basis of scale analysis. In DWT wavelet processing, is the approximate wavelet decomposition term, is the detailed wavelet decomposition terms in the x-direction and y-direction, and is the detailed wavelet decomposition terms in the diagonal direction. The subscripts of the terms represent the order of wavelet decomposition. During inverse DWT, only four–six order terms were included in this investigation. Consequently, the outline of the representative structures was distinguished by binarization with threshold processing, i.e., logical matrix 1 (). Application of the deconvolution method shrank the outline and further enhanced the spatial resolution, i.e., logical matrix 3 (). The image was then processed with repeated binarization, threshold analysis, and Gaussian interpolation before finalization. The overall processed image can be obtained by .

Figure 2

Schematic of DWDC method.

2.2.1. Expansion of Image by Gaussian Interpolation

To increase the resolution of the image, we first reduced the dot pitch. 3D Gaussian interpolation [37,38] was applied to expand the size of the image with reduced pitch. The original 20 images of 512 × 512 pixels were expanded twice to 80 expansion images () of 2048 × 2048 pixels. The Gaussian function is: where and are the lateral and axial Gaussian radius, respectively [39,40], is numerical aperture of lens, and is the wavelength of the excitation beam. After the expansion, the dot pitch was reduced to 63 nm.

2.2.2. Extraction of Microtubule Structures by DWT

In this investigation, we applied DWT to extract microtubule structures [41] from the expanded image matrix. Relative to the Fourier-transform-based image analysis, which can only filter images globally, DWT can extract different image structures based on their local scale and intensity distributions. This is more suitable for processing cellular microtubule structures, which are highly local and anisotropic. When performing wavelet decomposition, the information corresponding to the scale function is usually filtered by a low-pass filter, and the information corresponding to the wavelet function is filtered by a high-pass filter, as shown in Figure 3. At the same time, the scale information obtained by the low-pass filter can be used as the generating function of the wavelet function and the scale function of the next stage. The information corresponding to the scale function represents the low-frequency component in the original signal, which represents the coarse information of the original signal; the information corresponding to the wavelet function represents the high-frequency component in the original signal, which represents the detailed information component of the original signal.

Figure 3

Wavelet decomposition process, where is a low-pass filter, and is a high-pass filter.

For a two-dimensional image matrix, wavelet decomposition is processed in three directions, i.e., horizontal, vertical, and diagonal directions, as: where is the approximate wavelet coefficient after nth-order discrete wavelet decomposition, and , , and are the detailed wavelet coefficients of level n in horizontal, vertical, and diagonal directions, respectively, as shown in Figure 4. is the approximate information. is the detailed information in the x-direction, is the detailed information in the y-direction, and is the detailed information in the diagonal direction.

Figure 4

Diagram of second-order wavelet decomposition on . The left side is the image matrix for decomposition, and the right side shows the decomposed matrices. is the size of the image matrix. 2048 in this investigation.

The values are obtained from the image by the n-order DWT. , , , and have the following formats: where . After DWT, Equation (1) can be expressed as:

The image information can be further divided into two parts in scale space, i.e., the scales related to the desired structures (denoted as ) and the scales related to undesired structures (denoted as ), which can also be considered as noise structures. Thus, where where and correspond to the lower and higher bounds of DWT orders of the desired structures. For instance, if the desired structures have a characteristic size between 16 and 64 pixels, we have and . According to DWT decomposition, the image structures within desired scale ranges (or frequency ranges) can be retained in each local position. The undesired image structures, e.g., noise (which normally has high-frequency components) and image distortions due to nonuniform illumination (which has low-frequency components), can be removed from the images without the requirement of knowing their detailed distribution.

2.2.3. Binarization of Image

When we obtained the desired structures extracted by DWT, logic processing was carried out. The threshold value () was selected according to the probability density distribution of image intensity. The binarization image was thus obtained as a logical matrix:

After the processing, the image noise was further inhibited and the outline of demanded structure was highlighted.

2.2.4. Resolution Improvement by Lucy–Richardson (LR) Deconvolution Method

The logic matrix obtained in this way was not detailed enough and the structural resolution was not, apparently, improved. We then used the Lucy–Richardson deconvolution method to further process the image. Instead of directly applying LR on the image after DWT analysis, in this investigation, we applied LR on the logic matrix () to restore the sketch of the filament structures. This approach can prevent the spread of high-intensity structures that affect the low-intensity structures and lead to spurious images or overprocessing.

LR method is developed on the basis of Bayesian theory [42], Poisson distribution, and maximum likelihood estimation. The overall expression formula of LR deconvolution method algorithm is as follows: where is the iteration of logic matrix and are intermediate results. In each iteration of optimization, a scale factor is applied to evaluate the effect of this processing, according to the image difference before and after the iteration, as:

Thus, the iteration of the image during LR deconvolution can be obtained as: where is the maximum number of iterations.

During the deconvolution, one of the most important prerequisites is the evaluation of the actual PSF of the optical system. During LR deconvolution on the image, the PSF of the confocal imaging system is a two-dimensional (2D) Gaussian function, which can be expressed as , where is the actual PSF for LR deconvolution and is an experience coefficient that should be determined by a numerical experiment. However, since the DWT is followed by a binarization to extract the structure contour, the PSF of the image is broadened with a larger width than that of the original system. Furthermore, in the iterative process of deconvolution, background noise can also be gradually amplified, resulting in spurious images or overprocessing. Therefore, optimizing the number of iterations (i.e., ) and threshold deviations (damping coefficient) is important in the application of LR deconvolution. The evaluation of the image after deconvolution is highly arbitrary, depending on the desired structures of the image. In many cases, the restoration quality of an image cannot be simply evaluated by peak signal-to-noise ratio (PSNR) or structural similarity index measure (SSIM). Therefore, in this investigation, the deconvolution parameters were also determined by numerical experiments.

2.2.5. Post-Processing

After LR algorithm, the processed logic matrix showed nonlinear characteristics. Thus, we carried out another binarization process on , with the threshold value . The secondary logic matrix can be obtained as:

Finally, the final result () of DWDC processing can be obtained by multiplying the expanded image by , to extract the microtubule structure of the 3T3 cell, i.e.,

3. Experimental Results

3.1. Ground Truth Verification

The DWDC method algorithm was first tested by ground truth images, which are artificially generated bundles of filament-like structures, as shown in Figure 5a. The ground truth structures had a typical width of 63 nm. To mimic the actual biological image captured by a confocal microscope, a Gaussian blur and a Gaussian noise background were applied to the ground truth image, as shown in Figure 5b. The radius of the Gaussian blur function was 270 nm. The Gaussian noise had an average value of 0 and a standard deviation (STD) of 2. Before the DWDC processing, the blurred and noised image had a PSNR of 19.5773 dB and an SSIM of 0.0887 (see Table 1 for details). The latter was especially poor. After the DWDC processing, both the PSNR and SSIM increased to 20.6765 dB and 0.7082, respectively. The processed image in Figure 5c reproduces the skeleton of the ground truth with high visual consistency.

Figure 5

Comparison between the (a) ground truth image, (b) the blurred and noised image, and (c) the DWDC-processed image. Here, Gaussian noise with an average value of 0 and a standard deviation of 2 was applied to the ground truth to generate (b).

Table 1

PSNR and SSIM of the ground truth images applied with Gaussian noise and after DWDC processing.

Gaussian Noise		PSNR before Processing (dB)	PSNR after Processing (dB)	SSIM before Processing	SSIM after Processing
Average	STD
0	2	19.5773	20.6765	0.0887	0.7082
20	4	15.9681	19.7943	0.0425	0.4916
40	6	12.8995	18.7648	0.0330	0.3908
40	10	12.1049	17.8139	0.0297	0.3251

We further increased the noise level by applying a Gaussian noise with an average value of 40 and an STD of 10, to test whether a highly noised image could still be improved by the DWDC method. The results are presented in Figure 6. The strong noise (Figure 6b) broke the structure of the ground truth (Figure 6a) and led to misleading distributions, e.g., larger width and more peaks, as can be clearly observed from the profiles in Figure 7b. The PSNR and SSIM of the blurred and noised image are only 12.1049 dB and 0.0297. However, after the DWDC processing, both the PSNR and the SSIM increased to 17.8139 dB and 0.3251, respectively. The processed image in Figure 6c also reproduces the skeleton of the ground truth with acceptable visual consistency. The improvement is very clear.

Figure 6

Comparison between the (a) ground truth image, (b) the blurred and noised image, and (c) the DWDC processed image. Here, Gaussian noise with an average value of 40 and a standard deviation of 10 was applied to the ground truth to generate (b).

Figure 7

Comparison between the profiles in the ground truth image, blurred and noised image, and DWDC processed images; (a) corresponds to Figure 5, and (b) corresponds to Figure 6.

Figure 7 shows the profiles of the image intensity along a horizontal line under both the low noise (Figure 5) and high noise (Figure 6). In both cases, the profiles after DWDC processing showed significantly smaller FWHM than the blurred and noised profiles. In the low-noise case, the FWHM in the processed profile was 89.1 nm, which was close to that of the ground truth. In the high-noise case, the FWHM in the processed profile was 119.3 nm, which also approached that of the ground truth and exhibited clear resolution improvement relative to the blurred and noised image. It should be noted that although most of the noise peaks were removed, it is too difficult to remove all the noise in a highly noisy image, and some peaks may remain.

In Table 1, we list four pairs of comparisons before and after DWDC processing, with the different noise levels. In these cases, the PSNR improved obviously with greater noise and SSIM improved by a factor of nearly 10.

3.2. Expansion of Image by Gaussian Interpolation

A direct comparison between the original 512 × 512 image (i.e., matrix) and the expanded 2048 × 2048 image () is presented in Figure 8. First, visually, the expanded image shows the same high consistency as the original image, as shown in Figure 8a,b. At the corresponding positions, the intensity distributions along the selected row and column both show high similarity (Figure 8c,d), which shows that the details of the original image were reserved.

Figure 8

Comparison between the original and expanded images. White scale bars represent 10 µm. (a) Original image; (b) expanded image; (c) distributions of fluorescent intensity (NFI) on the corresponding horizontal and vertical positions before and after expansion; (d) zoom-in of (c).

3.3. Extraction of Microtubule Structures by DWT

In this investigation, we used Coiflet 3 wavelet function to decompose the image up to the sixth order, where 3 represents the subtype of the Coiflet wavelet. Because the microtubule structure in the image had a width of 20 to 60 pixels, in the extraction, we only kept 4–6 order components, i.e.,

In contrast to the expanded image (Figure 8b), the reconstructed image shown in Figure 9a clearly reserves the filament-like microtubule structure. The undesired components (Figure 9b), which make the image blurry and noisy, were successfully removed. After DWT analysis, a binarization process was applied, according to the probability distribution of , to extract the sketch of the desired structures. The probability distribution of is plotted in Figure 9c. Following the numerical experiments, we only retained the top 15% of the image intensity. Thus, the corresponding threshold value related to Figure 9c was estimated to be ; the sketch of the microtubule can be clearly observed in Figure 9d.

Figure 9

Extraction of microtubule structures by DWT. (a) Extracted components after DWT analysis. (b) Undesired components that should be abandoned. (c) Probability density function (PDF) of . (d) Logical matrix ( ) of the structure after DWT analysis. The white scale bars represent 10 µm.

3.4. Resolution Improvement by LR Deconvolution

At this stage, the sketch of the structure was still wide, and the spatial resolution of the structures was below our expectations (Figure 9d). Subsequently, we used the LR deconvolution algorithm to further process the logical matrix. During the numerical experiments, we used different , maximum iteration , and damping coefficient to optimize the outcome (i.e., narrow and consistent structural features). The representative results are listed in Figure 10a–d. Notably, all the images in Figure 10 were binarized after LR deconvolution.

Figure 10

Different empirical coefficients corresponding to the deconvolution results. White scale bars represent 10 µm. (a) ; (b) ; (c) ; (d) . (e) The comparison between pre-processing and post-processing of the logic matrix used to extract the microtubule structure in the original image. Here, is applied. The white scale bars represent 10 µm.

When , the applied Gaussian radium for deconvolution was only half of the theoretical value . The continuous structures of the logical matrix became fragmented and hollow (Figure 10a). When was increased to 2.5, significant shrinkage of the logical matrix structures was achieved with at the expense of continuous network structures (Figure 10b). The optimal outcome was obtained when and the damping coefficient were 10 and 0.01, respectively (Figure 10e). It should be noted that before applying the DWDC method, the optimal ξ must be calculated, since a mismatching ξ can lead to invalid or even over-processed images. For different optical systems, the PSF and the optimal ξ are different. For a different optical system, it may take a few hours to locate the optimal ξ.

3.5. Image after Processing

By extracting the structure of the logical matrix, the microtubule structure of 3T3 cells was obtained (Equation (12) and Figure 11a). The final image rendered the clear filament-like mesh structures of the microtubule in the 3T3 cells, with good contrast and low noise, as shown in Figure 11b. A clearer comparison between the original images and the processed images can be found in the 3D microtubule structures of the 3T3 cells in the Supplementary videos.

Figure 11

Comparison of original image and image processed after DWDC method. (a) Direct comparison before and after the DWDC method. The white scale bar represents 10 µm. (b) Local processing results of the DWDC method in the box of the dashed line of (a) for comparison with Figure 1b (see more details for 3D structure of Figure 11b in Supplementary video S2).

In Figure 12, we compare the distributions of the fluorescence intensity at the same positions of the original and processed images, which show a fifteen-fold improvement in spatial resolution, from 1.94 to 123.7 nm, as evaluated by the full width at half maximum (FWHM) of the structure. The image after DWDC processing achieved a super-resolution level, i.e., beyond the optical diffraction limit.

Figure 12

Comparison of FWHM between original image and image processing after DWDC method in the box of the solid line in Figure 11 (a). and are the image intensities of the 3T3 cell image; (b) is the zoom-in of (a).

Notably, the PSNR, a commonly used criterion for evaluating the noise level of images before and after processing, was −48.8. The SSIM, which is another widely used parameter, was only 0.0155. This implies that PSNR and SSIM, as the common judging standards, may not be the gold standard with which to evaluate feature extraction algorithms in biomedical and biological applications if drawing a comparison between original and processed images directly.

The difference before and after the DWDC processing from the 3D reconstruction can be seen more clearly. In Figure 13a, which is reconstructed from the original images, only the surfaces of the cells can be distinguished, not the microtubule structures. After processing, the previously unclear images of the microtubule network structure (Figure 1b,c) became considerably more distinguishable (Figure 13b), and more biological information was revealed. As an example, our results demonstrate that two cells (i.e., the green- and purple-colored cells) formed a cell–cell connection. When the third cell (yellow) passed through the gap between these two cells (Figure 13c), this remodeled its microtubule network and indicated that mechanical force was induced by the cell–cell collision [43]. Since it is widely accepted that propagating mechanical cues during the collective movement of population cells can activate mechano-signaling and regulate cellular behavior [44], in which the remodeling of cytoskeleton networks plays important roles, our approach shows potential for deciphering dynamic cytoskeleton network reorganization and remodeling at the single-molecule level, even with conventional high-resolution imaging techniques.

Figure 13

Schematic of cell collision and crossing through the 3D reconstruction of the images. (a) 3D reconstruction from the original images, (b) 3D reconstruction from the processed images, and (c) schematic of the relationship between the three cells.

4. Conclusions

In this investigation, we introduce the use of the DWDC method, which was developed based on the discrete wavelet transform and Lucy–Richardson algorithm, to extract the microtubule structures of 3T3 cells from confocal images. Using the DWDC method, a sequence of image processing steps was applied, including Gaussian interpolation, discrete wavelet transform, Lucy–Richardson deconvolution, binarization, and probability density analysis. The skeletons of the filament-like microtubule structures were distinguished with significantly reduced width. Finally, the microtubule structures can be extracted with much higher spatial resolution. The microtubule structure in the original image, which had a FWHM of up to 1.94 µm, was reduced to 123.7 nm after processing with the DWDC method. The improvement in the structural resolution was around fifteen-fold. In the numerical experiments, the PSNR of the original images was enhanced by up to 5.7 dB and the SSIM was improved by a factor of 10. Compared with the single use of discrete wavelet transform or Lucy–Richardson algorithm for image processing, the composite image processing method can effectively remove noise, improve the SNR, and increase the resolution of the image to a super-resolution level simultaneously. This investigation shows a new and effective approach to improving image resolution and SNR. It can be applied to wide-field and confocal microscopes, which are restricted by the optical diffraction limit, as well as to super-resolution microscopes, to further improve their imaging performance.