Ultrasound Speckle Reduction Based on Histogram Curve Matching and Region Growing

: The quality of ultrasound scanning images is usually damaged by speckle noise. This paper proposes a method based on local statistics extracted from a histogram to reduce ultrasound speckle through a region growing algorithm. Unlike single statistical moment-based speckle reduction algorithms, this method adaptively smooths the speckle regions while preserving the margin and tissue structure to achieve high detectability. The criterion of a speckle region is defined by the similarity value obtained by matching the histogram of the current processing window and the reference window derived from the speckle region in advance. Then, according to the similarity value and tissue characteristics, the entire image is divided into several levels of speckle-content regions, and adaptive smoothing is performed based on these classification characteristics and the corresponding window size determined by the proposed region growing technique. Tests conducted from phantoms and in vivo images have shown very promising results after a quantitative and qualitative comparison with existing work.

Speckle reduction techniques include the compounding methods of combining images of different frequency contents [Wang, Hong, Wu et al. (2019)] of different spatial views [Wilhjelm, Jensen, Jespersen et al. (2004)] or through temporal filtering in the time domain ]. Speckle can also be reduced by adaptively filtering images according to the information on the speckle formation to smooth out speckle while preserving the structure [Chen, Yin, Flynn et al. (2003); Dutt and Greenleaf (1996)]. As another nonlinear technique for simultaneously performing contrast enhancement and noise reduction, the anisotropic diffusion method [Li and Zhang (2017); Wang and Hiller (2000); Yu and Acton (2002)] can model image filtering as controllable heat flow. It is well known that an echo signal with sufficient speckle detected by the envelope (i.e., a large number of randomly distributed scatterings with a smaller interval than the ultrasound wavelength) shows a Rayleigh distribution. The deviation from this special scattering condition has been previously modeled by the Rice distribution (used to explain the "structure" or nonrandom coherent components) and the K distribution (used to explain the deviations from the nonuniformity in the scatter distribution). Many researchers have observed statistical information to develop adaptive filtering to rule out speckle noise. Modern ultrasound imaging systems always apply logarithmic compression to convert the large dynamic range of input echo signals to match the lower dynamic range of the display device. This nonlinear processing will change the statistics of the input envelope signals and make these statistical moments unreliable in terms of the large variance in a small-size window. Dutt et al. [Dutt and Greenleaf (1996)] analyzed the effect of nonlinear processing by theoretically studying the distribution functions and used the statistics of logarithmic compressed echo images from the Kdistribution model as the echo envelope to derive a variance-dependent parameter that can be used to quantify the extent of speckle formation. Instead of using the above statistical moments, we use the histogram shape to measure the similarity between the probability distribution of the running window and that from a user-selected reference window within the speckle region. The entire histogram shape should contain more information than single statistical moments. Furthermore, the whole histogram covers the changes in the distribution caused by system-dependent processing. In combination with other local image features, the estimated similarity value can be used to determine the region for adaptive averaging. Liu et al. [Liu, Czenszak and Kim (1998)] proposed using the histogram shape as the criterion to differentiate speckle regions, and this approach was extended in Li et al. [Li and Zhang (2017)]. In this paper, we propose a new method to define the processing window and calculate the similarity. Additionally, a fast algorithm is proposed to select regions for adaptive smoothing. Consequently, the basic ideas of our proposed algorithm are described as follows: (1) The local amplitude histograms are used to distinguish the speckle and non-speckle areas instead of the local statistical moments.
(2) A region growing method is adopted to select the optimal region of a local processing area rather than selecting a window with a constant size. This optimal region might be irregular and obtained from processing based on the similarity factor and local tissue features.
(3) Pixels with the same range of local similarity values are collected with the proposed region growing method and categorized into different tissue types. This process can avoid premature blockage caused by intruders from different tissue types, which may result in small smoothing windows that are too small, and can prevent structure blurring from an excessively large processing area. These difficulties always arise from a processing window with a fixed size.
(4) Once the final processing area is determined, simple arithmetical mean filtering is used to compute the output pixel from its adjacent pixels with the same tissue type. This approach may provide the maximum tissue contrast but will not blur tissue borders. In the following sections, we first introduce the theory of using histogram curve matching to detect speckle and how to define the similarity function by extracting the statistics embedded in the histogram shape (Section 2.1). We then present the tissue characterization and region growing methods for speckle suppression in Sections 2.2 and 2.3. In Section 3.1, we try to simulate echo ultrasound signals to determine whether the local statistical moments are reliable for speckle detection. The proposed new speckle reduction technique is also applied to tissue-mimicking phantom images and in vivo ultrasound images for algorithm verification. We compare the results with other speckle reduction methods in Section 3.2.

Speckle detection based on histogram curve matching
From a given reference window located in a speckle area, k bins and n samples can be used to calculate the amplitude histogram. The histogram of the reference speckle region has a characteristic shape. Other histograms with different shapes are considered to deviate from fully developed speckle. Mathematically, the optimal number of bins to represent a Gaussian distribution with a 95% confidence level from n samples is [Li and Zhang (2017) The number of bins could be 8 for 49 n = , i.e., from a local 7×7 window. This window size should be larger than the speckle size on the image, which can be calculated from the imaging system parameters: where lateral S and axial S are the averaging speckle sizes in the lateral and axial directions, Dl and Da are the lateral and axial envelope sampling distances, respectively, fs is the axial envelope sampling frequency, c f is the center frequency, D is the aperture size of a linear transducer in the lateral direction, 0 z is the distance from the transducer to the focal zone, and BW is the bandwidth of the envelope spectrum. Eq. (1) could be used to approximately determine the number of bins of more practically distributed bins, such as the double exponential from the compressed images. Moreover, the number of bins derived from Eq. (1) is assumed for non-empty bins from which we equally group the whole range of intensity (e.g., 255 for an 8-bit image), namely, 32 bins to have approximately 8 nonzero bins to represent a discrete histogram in practice. A good histogram representation can help users find an "acceptable" reference histogram of speckle statistics. In this article, we define a "quality box" bounded by 9 horizontal bins and 15 vertical pixels ( Fig. 1) for the reference histogram in a 7×7 window. These numbers can be estimated theoretically from histograms in the speckle area, which means that the reference histogram that the user selected could be verified from the quality box. If the processing window contains some resolved structures, then its histogram curve always shows a bias with respect to the reference window. For example, the window is narrow for echo-free (or signal saturated) regions, wider for specular reflectivity areas, or even multimodal for the edges. These windows will result in a "lower similarity" with a reference histogram of fully developed speckle. Similarly, the histograms from different comparison objects are also different. For instance, the width of the histogram can provide insights into the type of scattering related to the tissue type. Therefore, a quantitative measurement, namely, the similarity value, is necessary for evaluating the similarity between the two histogram shapes. The vectors ref and tm are defined as the histogram from the reference speckle area and the current processing window with the size of M×N, respectively. The difference between these two histogram shapes h can be defined from the sum-absolute-difference (SAD) of these two vectors starting from their individual centroid ref C and tm C to the two ends of both vectors. Mathematically, there is where Range L and Range R are the maximum values of the "left size" and the "right size" of ref and tm, respectively: In Eqs. (4) and (5), the operators left size, LR, and right size RR are: Then, the similarity value is defined as: Note that the denominator 2 M N × × is used for normalization in Eq.
(3), and the procedures from Eqs. (4)-(8) are used to make the similarity value more independent of the gain adjustments that are implemented in modern ultrasound scanners. The similarity value S shows statistically how similar the current local window is to the referenced speckle characteristics. A large S indicates that the processing histogram is similar to the reference histogram; otherwise, the processing window is not speckle noise and may include more resolvable structures.

Tissue characterization based on similarity values
An ultrasound B-mode image usually covers diverse tissue structures, including edges, cysts, lesions and so on. Lesions, for example, from the diffuse scattering of tissues may be revealed in the speckle statistics and a change in the mean grey level. Moreover, if a local scattering density is not related to fully developed speckle, then this densitydependent speckle characteristic can be detected from the similarity value of Eq. (9), i.e., from the curve matching of gray level histograms. Assuming that the local statistics of pixels in the lesion area are continuous in terms of the change in its histogram shape, we may segment tissue types by grouping the derived similarity into the characterization function, i Type : where , x y I is the pixel value at ( ) , x y on the processing image and Si is a predefined threshold for the i-th category of the tissue type. Since i S lies in the range of (0, 1], theoretically, it is a probe-dependent parameter that can be determined from tissuemimicking phantom images. For system optimization, however, we can conduct in vivo experiments to collect pixels selected by the user and their corresponding local histograms to determine those thresholds for different exam types.

Speckle reduction using region growing
The first purpose of the proposed speckle reduction algorithm is to categorize each pixel into one of the tissue types defined in Eq. (10). The collection of pixels belonging to the same tissue type in a concatenated region is called a homogeneous area. Therefore, each pixel defines its own homogeneous area, and the intended speckle reduction will be performed by smoothing each pixel's homogeneous area. Obviously, this homogeneous area is connected but could be very irregular. From the viewpoint of system implementation, the detection of an irregular shape requiring smoothing is time-consuming and difficult to execute in real time. Therefore, we propose an adjustable rectangular window that approximates the homogeneous region by reducing or increasing the size of the window, which depends on a function of the number of pixels with the same tissue type as the seed pixel (i.e., the current processing pixel). Note that pixels with the same tissue type are over the local window in terms of the seed pixel under the assumption of a continuous change in similarity values. We now summarize our speckle reduction algorithm as follows and give detailed descriptions in Section 2.4. Algorithm: Speckle Reduction Based on Histogram Curve Matching and Region Growing (1) Calculate the similarity value of Eq. (9) for every pixel from a local window of M N × and assign each pixel belonging to a tissue type i defined in Eq. (10) where T N is the total number of tissue types for one exam; (2) Define a set of smoothing windows of where NS represents the total number of windows to be smoothed; (3) Starting from the seed pixel ( ) 0 , I x y with a tissue type of 0 Step 3 for a new seed pixel until the entire image has been processed.

System implementation
In general, the clinical application of ultrasound is performed using different types of probes. For example, we can set Step 1 of the algorithm to perform a more detailed speckle classification in the abdominal examination from a 3.5 MHz curved array, and we can set 4 T N = for the shallow depth small parts scanned from a 7.5 MHz linear array. A setting of 7 T N = means that up to 7 different sizes of windows can be used for smoothing. This setting is a tradeoff between the level of smoothing and processing speed. In Step 2 of the proposed algorithm, we can define a set of windows with the size of 3 3, 5 5, 7 7, 9 9, 11 11, 13 13, 15 15 for later smoothing. Note that the number of tissue type NT is always larger than or equal to the number of smoothing windows NS. The calculation of the histogram of the tissue type j H in Step 3 of the proposed algorithm is for speeding up the process of counting the number of pixels in the processing window with the same tissue type as the seed pixel. During region shrinkage or growth, the counting of the seed tissue type can be obtained quickly by keeping the histogram vector to avoid deleting or adding the outermost pixels of the processing window.
The constant η used in Steps 4 and 5 of the proposed algorithm is to determine if the number of pixels that have the same tissue type as the seed pixel is dominant in the processing window. For instance, means that the number of pixels whose assigned tissue type is the same as the seed pixel is less than 70% of the total j j M N pixels of the processing window.
The constant η used in Steps 4 and 5 of the proposed algorithm is used to decide if the number of pixels that have the same tissue type as the seed pixel is dominant in the processing window. For instance, means that the number of pixels whose assigned tissue type is the same as the seed pixel is less than 70% of the total j j M N pixels. Moreover, the constant η could be extended according to the type of tissue. For example, we can set a larger η for the low level of tissue type (i.e., smaller similarity values) so that a smaller size of smoothing window is used to preserve the non-speckle tissue, such as the edges. Similarly, setting a smaller η for a higher level of tissue type will make region growing easier to occur, which is equivalent to having more smoothing in the speckle area.

Testing and discussion
Statistical moments are used in almost all the statistics-based image processing methods for speckle reduction, including the ratio of the mean to the standard deviation, the variance, and the ratio of the mean to the variance. These moments basically serve as criteria to check if the processed pixel could be classified as "speckle". In this section, it is verified that it is not reliable to use statistical moments obtained from a local window of a limited size as the criterion for speckle detection. Moreover, a well-known rule that "the statistics of pre-compressed envelope signals within a fully developed speckle area has a Rayleigh distribution and its signal-to-noise ratio (SNR) is defined as the ratio of the mean to the standard deviation of 1.91" may not be valid from the viewpoint of the discrete domain of the image processing. This rule means that an estimated SNR from the processing window with a limited size could be close to 1.91, even if the statistics of the signals in this local window are far away from the Rayleigh distribution. The histogram shape offers more information than the histogram statistics, can be used for speckle detection and performs well in actual situations.

Analysis of statistical moments and the histogram curve
Ultrasound images come from the envelope (or the amplitude) of the raw radio-frequency (RF) data and has no phase information. Without loss of generality, we now directly test the statistical moments from its random number generator (i.e., uncorrelated data) instead of using a convolution of the pulses and random scatters (correlated data) to generate the windowed pixels. To obtain the random deviation from the given probability distribution function (pdf), we can use the transformation method [Gonzalez and Woods (2002)]: Given a pdf ( ) By using the Rayleigh random number generator of Eq. (12), we can easily generate the datasets with different combinations of pdfs. Fig. 2 shows the histogram of a Rayleigh random generator where we generate 100 datasets with 1,000 numbers in each dataset. The histogram is the average of 100 histograms in Fig. 2, so that the statistics of the estimators can be ensured accurately. Fig. 3 records the statistical moments in 100 datasets, and the corresponding SNR is 1.9125, which is close to the theoretical value of a Rayleigh distribution of 1.91. Note that the calculated kurtosis has a large variance in a set of 1,000 numbers; thus, it could be an unreliable parameter in general windowed processing (due to the limited number of pixels). The histogram and statistical moments of a uniform distribution are shown in Figs. 2-5, where the computed kurtosis of -1.2 matches the theoretical value. From Figs. 2-5, it seems that the computed SNR from very different statistical distributions can show comparable values. This phenomenon may indicate that it is unreliable to use computed statistical moments to identify the inherent statistics (e.g., the distribution). Any misidentification will become serious when computed statistical moments come from pixels of limited size in an imaging window. In our application, we usually handle a small imaging window (for instance, 7×7), where the data length is 8 bits. We then compute the histogram from 49 numbers (Fig. 5), where the pixel amplitude has been grouped into 32 bins, as discussed in Section 2.1. The statistical moments from this small dataset show a large variance, as shown in Fig. 7. Similar results can be observed in the uniform distributions illustrated in Figs. 8 and 9. From Figs. 2-9, we can see that it is not reliable to use statistical moments as the criterion for image segmentation; we then verify our argument of better segmentation by using the proposed similarity values of histogram shapes. We can set up the data from a combination of the Rayleigh and uniform distributions by taking a certain percentage of our samples from one distribution and the remaining from the other distribution, i.e., 90%, 80%, 70%, 60%, and 50% from the Rayleigh distribution, and 10%, 20%, 30%, 40%, and 50% from the uniform distribution. The motivation for such a setup is the area of an image where the speckle borders significant electronic noise. Fig. 10 shows the histogram of the Rayleigh and uniform distributions since we tested 20 datasets, each with 49 numbers. Obviously, if more pixels come from the uniform distribution, then the histogram shape becomes flatter. A reference histogram can be obtained by averaging 20 histogram curves from the pixels that originate from the 100% Rayleigh distribution, as shown with the red curve (i.e., the bottom one) in Fig. 10. Similarity values of other histogram curves with different combinations of Rayleigh and uniform distributions are shown in Fig. 11, where 0.9269 is the mean of 20 similarity values relative to the reference histogram in the case of an 80% Rayleigh distribution. Moreover, we can compute the SNR from each 49 pixel set, where 1.63 is the mean of 20 SNRs from the 80% Rayleigh distribution. The standard deviation in Tab. 1 is 0.35, and the similarity value is recorded in Fig. 11. As shown in Tab. 1, the bias measurement is defined as the relative error of either the SNR or similarity value of a certain level of the Rayleigh distribution compared to the reference of the 100% Rayleigh distribution. It can be seen clearly from Tab. 1 that using the similarity value (i.e., the information from the histogram shape) as the criterion to detect how a certain level of the Rayleigh distribution deviates from the fully developed speckle (in our application, the ideal Rayleigh distribution) is much better than using the statistical moment (i.e., the SNR). This finding has been confirmed from the very high bias in the similarity compared with that in the SNR when the governing distribution increasingly deviates from the ideal Rayleigh distribution.

Experimental results from ultrasound signals
We used Saset Healthcare's iMago c21, a commercial digital ultrasound scanner designed and built by Saset (Chengdu) Inc., China for data acquisition. Fig. 12(a) shows the unprocessed phantom images from a 5 MHz linear probe scanned at a depth of 8 cm. Fig.  12(b) is the speckle reduction result from Dutt et al.'s method [Dutt and Greenleaf (1996)] based on a statistical moment (i.e., the variance) with a 7 7 × processing window for logcompressed ultrasound images. The result of our proposed algorithm is shown in Fig.   12(c), where the reference histogram is derived from a 7 7 × window selected manually in the speckle area from the original image. It is clear that, from Fig. 12, our result shows much better contrast resolution with largely smoothed speckle and preserved tissue structures (i.e., the point targets and a cyst). We then acquired in vivo images to verify our techniques. Fig. 13(a) is the human neck image from a 7.5 MHz linear array scanning at a depth of 4 cm. Fig. 13(b) shows the result from Dutt et al.'s method [Dutt and Greenleaf (1996)], and Fig. 13(c) demonstrates the speckle reduced image via our proposed method. In addition to linear scanning, we also tested a human liver image from a 3.5 MHz curved array ( Fig. 14(a)) and speckle reduced image from Dutt et al.'s method and our proposed method (Figs. 14(b) and 14(c)), respectively. The speckle reduced images from our method visually show better contrast resolution where the detailed tissue structure features were basically preserved without missing useful clinical information. Furthermore, to analyze the performance of the proposed speckle reduction algorithm in terms of the contrast improvement, we calculate the contrast-to-noise ratio (CNR) as [Chen, Yin, Flynn et al. (2003)]: Here, µc denotes the mean of the "object", µs is the average value of the background speckle that encloses the object, and σ 2 c and σ 2 s are variances of the object and the speckle areas, respectively. In general, the image area with a larger CNR shows better contrast resolution. Tab. 2 lists the CNR measurements where each value is an arithmetical mean of CNRs from the five selected areas with identified objects and speckle pixels large enough for comparison. It can clearly be seen from Tab. 2 that our proposed method always displays the largest CNR among the three sets of test images, which is consistent with the visual grayscale results presented in Figs. 12-14. This finding means that our techniques are promising for reducing ultrasound speckle.

Conclusions
We have proposed an ultrasound speckle reduction method based on local histogram curve matching and region growing filtering in the hope of adaptively reducing speckle noise while keeping the structure region unchanged. We have applied a SAD-based histogram matching technique of extracting the embedded statistical information from an image processing window with a limited size to differentiate ultrasound speckle from tissue structure. This scheme will output a similarity value from histogram matching, which can be used to adaptively classify the current region into a certain level of tissue type. We then smooth this processing area with the optimal region growing according to the classified tissue type. Since smoothing always occurs within the same tissue type, our proposed method not only offers a significant average value to reduce the local variance to improve the SNR but also preserves structure differentiation, which always shows different tissue types. This property has been verified by the CNR of the improved testing images.
In this article, we tested ultrasound phantom images and in vivo images. The results have shown that better contrast resolution can be achieved with our approach than with other speckle reduction methods, while the structure can also be preserved without losing useful clinical information. Currently, we are studying ways to speed up processing both on DSP  and GPU [Cook (2012)] platforms. More in vivo images are needed to fine-tune algorithms and parameters optimized specifically for the setting of feature patterns.

Conflicts of Interest:
The authors declare that they have no conflicts of interest in this research report.