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.


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 quasiisotropic features, e.g., round-shaped and nanometer-sized exosomes, a deconvolutionbased 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.

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.

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.
Micromachines 2022, 13, x 3 of 16 by the noise distribution. It is necessary to remove as much noise as possible before analyzing biological activities. 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).

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 • , .

Methods and Process
An optical image is a convolution of an object with the point spread function (PSF) of an optical system [36]. If M is the matrix of the image, where P is the PSF of the optical system, S is the light distribution according to the object and N 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, LL n is the approximate wavelet decomposition term, LH n is the detailed wavelet decomposition terms in the x-direction and y-direction, and HH n 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 (M e,dL ). Application of the deconvolution method shrank the outline and further enhanced the spatial resolution, i.e., logical matrix 3 (M e,LRL ). The image was then processed with repeated binarization, threshold analysis, and Gaussian interpolation before finalization. The overall processed image M F can be obtained by M e ·M e,LRL . 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 ( )

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 (M e ) of 2048 × 2048 pixels. The Gaussian function is: where r ⊥ = (0.61λ e ) N A and r = (4nλ e ) 2N A 2 are the lateral and axial Gaussian radius, respectively [39,40], N A is numerical aperture of lens, and λ e is the wavelength of the excitation beam. After the expansion, the dot pitch was reduced to 63 nm.

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.

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. Th original 20 images of 512 × 512 pixels were expanded twice to 80 expansion images ( of 2048 × 2048 pixels. The Gaussian function is: where = (0.61 ) ⁄ and ∥ = (4 ) ⁄ (2 ) 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.

Extraction of Microtubule Structures by DWT
In this investigation, we applied DWT to extract microtubule structures [41] from th expanded image matrix. Relative to the Fourier-transform-based image analysis, whic can only filter images globally, DWT can extract different image structures based on thei local scale and intensity distributions. This is more suitable for processing cellular micro tubule structures, which are highly local and anisotropic. When performing wavelet de composition, the information corresponding to the scale function is usually filtered by low-pass filter, and the information corresponding to the wavelet function is filtered by high-pass filter, as shown in Figure 3. At the same time, the scale information obtained b the low-pass filter can be used as the generating function of the wavelet function and th scale function of the next stage. The information corresponding to the scale function rep resents the low-frequency component in the original signal, which represents the coars information of the original signal; the information corresponding to the wavelet function represents the high-frequency component in the original signal, which represents the de tailed information component of the original signal.  For a two-dimensional image matrix, wavelet decomposition is processed in three directions, i.e., horizontal, vertical, and diagonal directions, as: where c n is the approximate wavelet coefficient after nth-order discrete wavelet decomposition, and d n,x , d n,y , and d n,D are the detailed wavelet coefficients of level n in horizontal, vertical, and diagonal directions, respectively, as shown in Figure 4. m n,LL is the approximate information. m n,LH is the detailed information in the x-direction, m n,HL is the detailed information in the y-direction, and m n,HH is the detailed information in the diagonal direction.

, ,
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. 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,  The values are obtained from the image M e by the n-order DWT. m n,LL , m n,LH , m n,HL , and m n,HH have the following formats: where m 0,LL (x, y) = M e (x, y). 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 M e,d (x, y)) and the scales related to undesired structures (denoted as M e,ud (x, y)), which can also be considered as noise structures. Thus, where M e,d (x, y) = ∑ n=n h n=n l d n,x m n,LH (x, y) + d n,y m n,HL (x, y) + d n,D m n,HH (x, y) where n l and n h 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 n l = 4 and n h = 6. 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 highfrequency 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.

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.

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 (M e,dL ) to restore the sketch of the filament structures. This approach can prevent the spread of highintensity 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 M k (x, y) is the kth iteration of logic matrix M e,dL and M mid (x, y) are intermediate results.
In each iteration of optimization, a scale factor f is applied to evaluate the effect of this processing, according to the image difference before and after the iteration, as: Thus, the kth iteration of the image during LR deconvolution can be obtained as: where kmax 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 P = P 0 exp −[(x − x c ) 2 + (y − y c ) 2 /2(∆r) 2 , where ∆r = ξr ⊥ 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., kmax) 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.

Post-Processing
After LR algorithm, the processed logic matrix showed nonlinear characteristics. Thus, we carried out another binarization process on M kmax , with the threshold value χ = 1. The secondary logic matrix can be obtained as: Finally, the final result (M F (x, y)) of DWDC processing can be obtained by multiplying the expanded image M e by M e,LRL , to extract the microtubule structure of the 3T3 cell, i.e.,

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.

Ground Truth Verification
The DWDC method algorithm was first tested by ground truth images, which artificially generated bundles of filament-like structures, as shown in Figure 5a. ground truth structures had a typical width of 63 nm. To mimic the actual biological im captured by a confocal microscope, a Gaussian blur and a Gaussian noise backgro were applied to the ground truth image, as shown in Figure 5b. The radius of the Gaus blur function was 270 nm. The Gaussian noise had an average value of 0 and a stand deviation (STD) of 2. Before the DWDC processing, the blurred and noised image ha PSNR of 19.5773 dB and an SSIM of 0.0887 (see Table 1 for details). The latter was e cially poor. After the DWDC processing, both the PSNR and SSIM increased to 20.6765 and 0.7082, respectively. The processed image in Figure 5c reproduces the skeleton of ground truth with high visual consistency.   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. broke the structure of the ground truth ( Figure 6a) and led to misleading distrib e.g., larger width and more peaks, as can be clearly observed from the profiles in 7b. The PSNR and SSIM of the blurred and noised image are only 12.1049 dB and However, after the DWDC processing, both the PSNR and the SSIM increased to dB and 0.3251, respectively. The processed image in Figure 6c also reproduces the s of the ground truth with acceptable visual consistency. The improvement is very  Figure 7 shows the profiles of the image intensity along a horizontal line und the low noise ( Figure 5) and high noise ( Figure 6). In both cases, the profiles after processing showed significantly smaller FWHM than the blurred and noised pro the low-noise case, the FWHM in the processed profile was 89.1 nm, which was that of the ground truth. In the high-noise case, the FWHM in the processed prof 119.3 nm, which also approached that of the ground truth and exhibited clear res improvement relative to the blurred and noised image. It should be noted that al most of the noise peaks were removed, it is too difficult to remove all the noise in a 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.

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 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.

Expansion of Image by Gaussian Interpolation
A direct comparison between the original 512 × 512 image (i.e., M matrix) and the expanded 2048 × 2048 image (M e ) 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.
A direct comparison between the original 512 × 512 image (i.e., matrix) and th expanded 2048 × 2048 image ( ) is presented in Figure 8. First, visually, the expande image shows the same high consistency as the original image, as shown in Figure 8a,b. A the corresponding positions, the intensity distributions along the selected row and co umn both show high similarity (Figure 8c,d), which shows that the details of the origina image were reserved.

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 M e,d , to extract the sketch of the desired structures. The probability distribution of M e,d 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 χ = 61; the sketch of the microtubule can be clearly observed in Figure 9d.

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 kmax, 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. ure 9a clearly reserves the filament-like microtubule structure. The undesired co (Figure 9b), which make the image blurry and noisy, were successfully remo DWT analysis, a binarization process was applied, according to the probabilit tion of , , to extract the sketch of the desired structures. The probability dist , is plotted in Figure 9c. Following the numerical experiments, we only re top 15% of the image intensity. Thus, the corresponding threshold value related 9c was estimated to be = 61; the sketch of the microtubule can be clearly o Figure 9d.

Resolution Improvement by LR Deconvolution
At this stage, the sketch of the structure was still wide, and the spatial re the structures was below our expectations ( Figure 9d). Subsequently, we used convolution algorithm to further process the logical matrix. During the numeri ments, we used different , maximum iteration , and damping coefficie When ξ = 0.5, the applied Gaussian radium for deconvolution ∆r was only half of the theoretical value r ⊥ . 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 kmax 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 ξ. 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 ξ.

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.
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 µm 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.
was obtained (Equation (12) and Figure 11a). The final image rendered the clear filamentlike 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.
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 μm 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.  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 purplecolored 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. Micromachines 2022, 13, x 13 of 16 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.

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 fifteenfold. 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.

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.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/mi13060824/s1. Supplementary video S1: three-dimensional reconstruction of the microtubule skeleton of 3T3 cells (before processing); supplementary video S2: three-dimensional reconstruction of the microtubule skeleton of 3T3 cells (after processing).