Speckle noise reduction in coherent imaging systems via hybrid median–mean filter

Abstract. Images recorded by coherent imaging systems, including laser-based photography, digital holography (DH), and digital holographic microscopy (DHM), are severely distorted by speckle noise. We present a single-shot image processing method to reduce the speckle noise, coined hybrid median–mean filter (HM2F), which is based on the average of conventional median-filtered images with different kernel sizes. The synergic combination of the median filter and mean approach provides a denoised image with reduced speckle contrast while the spatial resolution is kept at up to 97% of the original value. The HM2F method is compared with the conventional median filter approach, the 3D block matching filter, the nonlocal means filter, the 2D windowed Fourier transform filter, and the Wiener filter using different speckle-distorted images to benchmark its performance. Based on the experimental results and the simplicity of the technique, HM2F is proposed as an effective denoising tool for reducing the speckle noise in laser-based photography, DH, and DHM.


Introduction
Speckle noise is inherited by all imaging techniques that utilize a coherent light as a source of illumination. The speckle is generated due to the illumination of a surface with roughness in dimensions comparable to the coherent wavelength of the light source. 1 Following Huygens' principle, 2 each point of the surface acts as a spherical wave scatter emitting light with random phase distribution among them. The coherent superposition of these spherical wavefronts creates a random intensity pattern of constructive and destructive interference, dark and shiny spots, to produce the speckle pattern. 1 The speckle noise is, therefore, a random pattern of dark and shiny spots that appears on the images in any coherent (e.g., laser-based) imaging system, 1 including laser-based photography, 3 digital holography (DH), 4 and digital holographic microscopy (DHM). 5 While laser-based projectors provide a wide color gamut for vivid, super bright, and high contrast images, their major limitation is the speckle noise caused by the lasers' coherent nature. Akram and Chen 3 reviewed the current state-of-the-art solutions to reduce speckle noise. These approaches involve the decorrelation of the light through the wavelength using a source with the fast tuning of lasing wavelength or an array of sources with different wavelengths. Another decorrelation approach can involve the use of a source with reduced spatial coherence (e.g., spatial decorrelation) or the use of an array of sources with different angles at one spot on the screen (e.g., angular decorrelation). The angular decorrelation can also be achieved changing the illumination angle at one spot or rotating a diffuser. Finally, light can be decorrelated using polarization with a source with fast changing polarization or splitting the light into two paths with sufficient path length difference and different polarizations. *Address all correspondence to Ana Doblas, adoblas@memphis.edu Although DH and DHM techniques have been successfully applied to many diverse fields, including life and material sciences, [5][6][7] the presence of speckle noise has hampered additional potential applications. Several approaches have been developed in DH and DHM to diminish the adverse speckle effects to widen the use of these imaging methods; however, these ruin the image contrast and reduce the spatial resolution. In general, the methods to reduce speckle noise can be grouped into two categories: physical approaches [8][9][10][11][12][13][14][15][16] and computational approaches. [17][18][19][20][21][22][23][24][25][26][27] The physical techniques are commonly based on two methodologies. One method consists of recording the same scene hologram Q times under different physical conditions. Because of the random behavior of the speckle noise, the speckle noise in each of the recorded holograms is uncorrelated with the others. Therefore, the average of the reconstructed images provides a numerical reconstruction with reduced speckle noise following a 1∕ ffiffiffiffi Q p law. 28 Such holograms can be acquired by changing the angle of illumination through a rotating mirror, 8 slightly rotating the object, 9 and lateral shifting the object 10 or the sensor. 11 Another reported physical technique consists of using different wavelengths to register multiple holograms. 12 In addition, other physical methods use some optical elements to reduce the light's spatial coherence partially. Examples of these elements include rotating diffusers, 13 spatial light modulators, 14,15 and holography diffusers. 16 Due to the number of holograms that need to be recorded to produce significate effects ranging from 20 to 100 images, the common disadvantage of all of these physical methods is the lengthy acquisition time, which limits the techniques to static samples. Computational techniques aim to reduce the speckle noise in dynamic applications. The main idea of these techniques is the use of filters applied in the spatial 17,21,23 or Fourier domain 24,25 to reduce the speckle noise. Other computational methods are based on the same idea as the physical approaches in which one averages Q partially uncorrelated speckle-distorted reconstructed images from a single hologram. These methods include the introduction of a spatial jittering in the Fresnel Kernel, 26 the use of a mask to generate Q sub-holograms, 27 and the sequential sampling of the discrete Fourier transform of the hologram. 19 A standard limitation of these computational approaches is the requirement of some parameters for improving the performance. This may be a limitation for inexperienced users who may not have prior knowledge of the technique. In addition, computational methods are subject to a trade-off between spatial resolution and speckle noise. The higher the reduction in the speckle noise is, the lower the spatial resolution in the denoised image is. This trade-off has been avoided by implementing approaches in which physical and computational methods are integrated. 19,29,30 In particular, Bianco et al. 30 proposed a framework that combines the acquisition of multiple digital holograms with optimized jointaction computational image-denoising methods. In that work, the authors demonstrated their technique using computer-generated holograms from a single-shot recording and validated it using dual-wavelength holographic recordings.
With the significant advancement of deep learning (DL) techniques over the last years, DL approaches have also been applied to restore sharp images from their degraded version in the presence of speckle noise. Among the different DL approaches, speckle-free images have been achieved using a convolutional neural network 31 and conditional generative adversarial network, 32 which were applied to underwater sonar images and laser-illuminated ex vivo porcine gastrointestinal tissues, respectively. Whereas traditional DL algorithms require the speckle-free image for training the model, Yin et al. 33 proposed a DL method that does not require prior knowledge of speckle-free object distribution. Even though DL algorithms are robust, their performance depends heavily on the number of training data and their quality.
The proposed method in this work synergistically combines two standard denoising methods in image processing: the median filter and the mean approach. The novelty of the presented approach, named the hybrid median-mean filter (HM2F), is the cascaded application of these techniques to the speckle-distorted images. In this work, we focus on the implementation of the HM2F method in speckle-distorted images from laser-based photography, DH, and DHM. Using the HM2F, the speckle noise is reduced up to 49% for laser-based photography, 72% for DH, and 14% for DHM, which corresponds to a speckle contrast of 0.51, 0.28, and 0.86, respectively. The resolution is kept at up to 97% of the original value, reducing the adding of blurring effects for all imaging modalities (e.g., laser-based photography, DH, and DHM).

Hybrid Median-Mean Filter
The HM2F method is an iterative method based on applying i times a median filter over the original noisy speckle-distorted image g. For each i iteration, the image after applying the median filter is saved in the variable h. The square kernel size of the median filter increases by (2i − 1) for each iteration. Because the median filter with kernel size ½1 × 1 provides an identical image to the noisy input image (g), the first iteration in the HM2F corresponds to i ¼ 2. The novelty of the HM2F to other methods reported in the literature is the combination of several median-filtered images (h) to provide a final denoising image (ĝ). We propose the average between the i'th median-filtered image (h) and the previous (i − 1)'th denoised image (ĝ), namely,ĝ. Therefore, in the first iteration (i ¼ 2), h ¼ medianðg; ½3 × 3Þ and g ¼ ðg þ hÞ∕2. Figure 1 shows the pseudocode of the HM2F. From Fig. 1, one can realize that the inputs of HM2F are the noisy images (g) and the maximum kernel size for the median filter (k). Figure 2 shows the process of the HM2F for the maximum kernel size of k ¼ 5 (maximum iteration, i ¼ 3). For k ¼ 5, the median filter is applied twice to the noisy image with kernel sizes of ½3 × 3 and ½5 × 5. For the first iteration, a median filter of kernel size ½3 × 3 is applied to the original noise image (g) to provide the denoised h image. The average of this result and the original image results is the new imageĝ. For the second iteration, a new denoised h image, obtained by applying the median filter with kernel ½5 × 5 to the noisy image, is averaged with the previously computed image. Figure 2 shows the operation of the median filter in two different regions, marked by the red-dashed boxes. Black values represent the original noisy data, and red values represent the result after applying the HM2F. A zero-padding operation is required to offset border problems for each median filter.
The HM2F can also be applied to RGB speckle-distorted images after splitting the color image into three-channel images, applying the HM2F to each channel image, and merging back these three-channel denoised images.

Experimental Validation of HM2F
The following section is devoted to showing the versatility of the proposed HM2F in laser-based photography, DH, and DHM. The HM2F was applied to the RGB noisy image of a four-color cube in laser-based photography, the reconstructed amplitude images of a dice and horse model in DH, and the reconstructed phase images of a star target and a complex biological sample in DHM. Figure 3 shows the optical configuration of the three optical systems. Figure 4 compares the performances of the HM2F method with the known denoising methods, including the 3D block matching (BM3D), 20 the nonlocal means (NLM) filter, 18,21 the Wiener filter, 22 and the conventional median filter (CMF). Note that for CMF we applied a median filter with a particular kernel size over the original image. For this comparison, we show the final RGB denoised image of a four-color cube utilized as a millimeter-sized object. In the laser-based RGB  photographic system [ Fig. 3(a)], the cube was directly illuminated by three lasers of wavelengths 671 (red), 532 (green), and 473 nm (blue). The cube was located at 86 cm from a commercial Canon photographic camera. More details of this laser-based photographic system are provided in Ref. 34. The kernel size for both the CMF and the proposed HM2F was 11 × 11. Using these results, we quantified and characterized the reduction in the speckle noise by measuring the speckle contrast (C) per channel as C ¼ σ i ∕I i , 28 where σ i and I i are the standard deviation and mean intensity, respectively, inside the yellow region marked in the original noisy image in Fig. 4(a). We report the average value of the three measured speckle contrast values for each denoising technique in the lower right corner of Fig. 4. These values were normalized to those of the original noisy image. Comparing the speckle contrast values, the reduction in the speckle contrast provided by the proposed HM2F method is almost the same as that provided by the CMF approach (C ¼ 0.51 for HM2F versus 0.50 for CMF). Nonetheless, for the same kernel size, the blurring effect of the conventional median-filtered RGB image is more significant than that of the HM2F image. A quantitative analysis of the blurring effects is provided by analyzing the first-order derivative of a step function (i.e., edge) for each method. The step function is defined by the mean profile along the vertical direction of the region marked by the green rectangle in Fig. 3(a). Because the first-order derivative of a step function is a Delta function, 35 we estimated the full width at half maximum (FWHM) of the first-order derivative of the step function. Because our step function is noisy, we fitted it using a smoothing spline with a parameter that was determined empirically for each method in such a way that the correlation coefficient between the noisy step function and the fitted one is higher than 0.99. Figure 4(g) shows the firstorder derivative of the edge for the BM3D, NLM, Wiener, CMF, and HM2F approaches. For each of these profiles, we quantified the FWHM value, reported in millimeters in the top right corner of each image method. The FWHM is equal to 0.28 mm for the original noisy image and the denoised images obtained after applying the BM3D, NLM, Wiener, and HM2F approaches. This means that the BM3D, NLM, Wiener, and HM2F methods do not introduce any blurring effect. By contrast, the CMF technique introduces some blurring with the FWHM being increased to 0.31 mm. Nonetheless, while the NLM and HM2F approaches are similar, the reduction of the speckle contrast is less for the BM3D (C ¼ 0.59 for BM3D versus 0.51 for HM2F) and the Wiener (C ¼ 0.58 for Wiener versus 0.51 for HM2F) methods. In conclusion, for the four-color cube, the HM2F shows the highest reduction the speckle contrast without adding blurry effects.

DH Results
The holograms for two different samples (dice and horse model) were recorded in the DH system [ Fig. 3(b)] operating in an off-axis Mach-Zehnder interferometer architecture. 36 The illumination source used in the DH system was a He-Ne laser of wavelength λ ¼ 632.8 nm. The holograms were recorded by a complementary metal-oxide-semiconductior (CMOS) camera (M × N ¼ 2592 × 2048 and Δx ¼ Δy ¼ 4.8 μm square pixels). In this experiment, the camera was located at a distance of z ¼ 70 cm from the sample, producing a speckle grain with a lateral (x − y) size of δ x ¼ λz∕MΔx ¼ 35 μm and δ y ¼ λz∕NΔy ¼ 45 μm. . We quantified and characterized the reduction in the speckle noise by plotting normalized speckle contrast versus the number of iterations for each kernel size (i.e., iteration) in Fig. 4(e). Because the speckle contrast is highly dependent on the object information, the contrast speckle was measured in 10 different square regions to provide an experimental error. The quantitative values shown in Fig. 5(e) show that the speckle contrast reduces rapidly for both the CMF and HM2F. Comparing the reconstructed amplitude images from the CMF and HM2F shows that, although the CMF reduces the speckle noise faster than the HM2F, the HM2F provides final denoising images with fewer blurring effects than that provided by the CMF. These results show that the HM2F presents a better trade-off between reducing the speckle contrast and the blurring effects. A quantitative analysis of the blurring effects is provided by estimating the FWHM of the first-order derivative of a step function defined by the profile along the transverse direction marked by the green line in Fig. 5(a). Figure 5(e) also compares the increase of the blurring effects versus the size of the kernel for the CMF and the HM2F. This panel illustrates how the speckle contrast and the response to the blurring can be controlled by selecting a different kernel size of the median filter. These curves show the trade-off between the speckle noise and resolution in both the CMF and HM2F methods. Despite this trade-off being present in both CMF and HM2F methods, the HM2F method shows a smaller blurring effect than the CMF method for any kernel size of the median filter, making it a more suitable method for mitigating the trade-off between speckle noise and blurry effects than the CMF approach.
To finalize the discussion for the DH system, Fig. 6 compares the performances of the proposed method with the known denoising methods [e.g., the BM3D filter, 20 the NLM filter, 18,21 the Wiener filter, 22 and the 2D windowed Fourier transform filter (WFT2F)]. 25 For this comparison, we used the reconstructed amplitude image for the two samples: dice (first row) and model horse (second row). Figure 6 shows that the HM2F performs equally well in high contrast images (i.e., dice results) and low contrast images (i.e., horse results). We measured the speckle contrast (C) inside the green region marked in the original images for the two samples. The speckle contrast values, which were normalized to those of the original noisy image, are reported in the lower right corner. As can be seen, the value of the speckle contrast varies from image to image. For example, for the dice image, the BM3D approach provides the denoised image with the lowest speckle contrast (i.e., C ¼ 0.20). The NLM and HM2F filters give simiar values on the speckle contrast. Table 1 reports the FWHM values of the first derivative of the step function defined by a transverse green line, depicted in Fig. 6(f). Based on the FWHM values, the Wiener, BM3D, and HM2F approaches produce the denoised amplitude images with the smallest blurring effect. Nonetheless, although the Wiener filter has a lower reduction of the resolution, the speckle contrast is even higher than that of the other methods. For this sample object, the speckle contrast of the denoised BM3D horse image is slightly lower than that of the denoised HM2F horse image (C ¼ 0.24 for the BM3D versus C ¼ 0.28 for the HM2F). Nonetheless, the HM2F method provides a final denoised amplitude image with improved contrast, as the histogram of Fig. 6(n) shows. Figures 6(m) and 6(n) show that the original noisy image and the denoised images for the dice and horse model present many pixels with low gray levels. The higher the amount of low gray levels is, the lower the mean value (I i ) is and, consequently, the higher the speckle contrast is because C is inversely proportional to I i . Nonetheless, the presence of low gray levels also affects the standard deviation value, resulting in low values of σ i . For each image in Fig. 6, we calculated the mean and standard deviation in 10 different zones of 50 × 50 pixels for the dice model and 20 × 20 pixels for the horse model. The average value of these parameters is reported below each panel in Fig. 6. Table 1 also reports the signal-noise ratio (SNR) of the original and denoising images for the dice and horse samples. The SNR value was estimated as 10 log 10 ðg max ∕stdðg b ÞÞ, where g max is  the maximum value of the amplitude image over the whole field of view and stdðg b Þ refers to the standard deviation of a region of the image in which there are no sample details. Because the SNR value depends on the region chosen to compute the standard deviation, we selected 10 different areas within the dice and horse field of views. Table 1 reports the average values of the SNR values for these two DH images.

DHM Results
The ability of the proposed method to reduce speckle contrast and mitigate the blurring effect was also demonstrated in quantitative phase imaging recorded in the DHM system. The optical setup of the off-axis DHM system follows a Mach-Zehnder architecture [ Fig. 3(c)], operating in a telecentric configuration 37,38 and image plane [e.g., the sensor was placed at the back focal plane of the tube lens (TL)]. Here, we validated the HM2F approach using two samples: a star target of the commercial quantitative phase target from Benchmark Technologies and a transverse section of the head of a Drosophila melanogaster fly. For the star target, the illumination source was a diode laser of wavelength 532 nm, and the CMOS camera had 1920 × 1080 square pixels with a 2.4 μm pixel size. The DHM image system consisted of an infinity-corrected 40×∕0.65 Olympus microscope objective (MO) and a TL of focal length 200 mm, resulting in an effective lateral magnification equal to 44.44×. On the other hand, the hologram of the Drosophila melanogaster fly was recorded using a HeNe laser of wavelength 633 nm and a charge-coupled device camera with 1024 × 1024 square pixels of 6.9 μm size. The imaging system was set up using an infinity-corrected 10 × ∕0.45 Nikon MO lens and a TL of focal length 200 mm. The star target allows us to measure the experimental resolution limit (RL) by estimating the minimum resolvable star pattern and quantifying how much the resolution was reduced for each approach. We defined the experimental RL as the diameter in which the contrast of the reconstructed phase star pattern was reduced by 10% from its reference value, which is the contrast value for the diameter equal to 150 μm. The results of the star target are shown in Fig. 7. For each method, the experimental RL is marked by a black-dashed circle. Whereas Fig. 7 Fig. 7(g)]. For the methods that use a median approach [Figs. 7(b) and 7(g)], the kernel size applied was ½5 × 5. Whereas the BM3D filter approach reduces the RL by a factor of 43×, the RL in the reconstructed phase image after applying the HM2F is not reduced. In addition to estimating the RL, we measured the speckle contrast (C) as a metric to calculate the reduction in the speckle noise. The parameter C was measured inside the yellow region marked in the inset of Fig. 7. The lower the value of C is, the higher the reduction in the speckle noise is. The reconstructed phase image obtained by the NLM filter [ Fig. 7(e)] is the most insensitive to the speckle noise (i.e., the smallest C value). However, the RL was highly diminished by 12%, reducing the ability to discriminate the finer details. By contrast, the HM2F provided the second smallest C value, resulting in the best method that minimized the trade-off between the reduction in the RL and the speckle noise.
Finally, Fig. 8 shows the reconstructed phase image of a section of the head of a D. melanogaster fly, a biological complex sample, provided by the BM3D, the NLM filter, the WFT2F, the Wiener filter, and the proposed HM2F. For comparison purposes, the noisy reconstructed phase image is shown in panel (a). For this sample, the BM3D and WFT2F approaches do not provide an optimal result. Nonetheless, the performance of the NLM, Wiener, and HM2F filters are similar. When an intense decorrelated noise corrupts the phase values, the reconstructed phase image has a high number of phase jumps in Fig. 8(g), marked by the white and black colors. We quantified the phase jumps in the original noisy image and denoised phase images inside the marked yellow region in Fig. 8(a). The phase jumps were reduced from 201 pixels in the original noisy image [ Fig. 8 Fig. 8(j)]. There is a reduction of 60% in the number of phase jumps for the HM2F, which shows that the proposed method provides a reconstructed phase image with reduced speckle noise.

Conclusions
In summary, we have demonstrated a single-shot denoising image method to reduce speckle noise. The HM2F is based on the synergetic combination of two well-known approaches in image processing: the median and mean filters. Our experimental results demonstrate that the HM2F reduces the speckle noise in color laser-based photography and both reconstructed amplitude and phase images from DH systems with a minimum addition of blurring effects. The performance of the HM2F approach was compared with that of the state-of-the-art methods in speckle denoising (e.g., BM3D, NLM, WFT2F, and Wiener). Based on our experimental results, the performance of the HM2F is more constant across the different types of images. For example, the WFT2F method does not provide satisfactory results for the dice and horse model [Figs. 6(d) and 6(j)]. By contrast, the NLM approach works for the dice image but fails for the horse model. The denoised Wiener image of the dice model presents more speckle noise (e.g., C > 1). Regarding the BM3D method, the denoised BM3D images for laser-based photography and DH present an adequate balance between the speckle reduction and the blurring effects. Nonetheless, this technique fails for quantitative phase imaging in DHM, reducing the RL by 43%. Another advantage of the HM2F is its simplicity in the required number of parameters and the computational processing time. While the HM2F and Wiener filter only require a single parameter (e.g., the kernel size) to provide a successful denoised image with reduced speckle noise, the WFT2F method requires the correct knowledge of the nine required parameters. Note that the maximum kernel size in the HM2F can be determined manually to obtain the best denoised HM2F image quickly and efficiently without requiring users' experience or knowledge. Regarding the processing times, the average processing time for the images of the dice and horse model and the star target is 1.4 s for the HM2F, whereas the times for the Wiener filter, BM3D, NLM, and WFT2F methods are around 0. 29,9,73, and 4390 times the processing times of the HM2F, respectively. The processing times are reported based on a 10 Pro Windows-based AMD Ryzen 5 8192 MB RAM laptop computer. Based on these results, the simplicity of the technique, and the processing time, HM2F is proposed as an effective denoising tool for reducing the speckle noise in laser-based photography, DH, and DHM. For increasing the use of our method to the community, the HM2F method was implemented as a script for Python and MATLAB and is available in a public repository on GitHub-https://oirl.github.io/Speckle-Hybrid-medianmean/. We predict that the HM2F could be applied to any image distorted by speckle noise. In future work, we will investigate the impact of this method on more interferometric systems, such as ultrasound and optical coherence tomography.
Jorge Garcia-Sucerquia received his PhD in physics from the Universidad Nacional de Colombia-Sede Medellin, Medellin, Colombia, in 2003. From 2004 to 2006, he held a postdoctoral position at Dalhousie University, Halifax, Canada, to work on digital lensless holographic microscopy. From 2012 to 2013, he was a visiting professor at the Universitat de València, Valencia, Spain. Currently, he is an associate professor at the School of Physics, Universidad Nacional de Colombia-Sede Medellin, Medellin, Colombia. His main research interests include DHM and their applications.
Ana Doblas is an assistant professor at the University of Memphis. She received her MS degree and PhD in physics from the Universitat de València, Valencia, Spain, in 2011 and 2015, respectively. Since 2019, she has been the principal investigator of the Optical Imaging Research Laboratory (OIRL). Her research interests include optical engineering, computational optics, and three-dimensional imaging, particularly designing novel microscopic imaging systems and their applications.