Hybrid Speckle Noise Reduction Method for Abdominal Circumference Segmentation of Fetal Ultrasound Images

Received Sep 17, 2017 Revised Mar 12, 2018 Accepted Mar 20, 2018 Fetal biometric size such as abdominal circumference (AC) is used to predict fetal weight or gestational age in ultrasound images. The automatic biometric measurement can improve efficiency in the ultrasonography examination workflow. The unclear boundaries of the abdomen image and the speckle noise presence are the challenges for the automated AC measurement techniques. The main problem to improve the accuracy of the automatic AC segmentation is how to remove noise while retaining the boundary features of objects. In this paper, we proposed a hybrid ultrasound image denoising framework which was a combination of spatial-based filtering method and multiresolution based method. In this technique, an ultrasound image was decomposed into subbands using wavelet transform. A thresholding technique and the anisotropic diffusion method were applied to the detail subbands, at the same time the bilateral filtering modified the approximation subband. The proposed denoising approach had the best performance in the edge preservation level and could improve the accuracy of the abdominal circumference segmentation. Keyword:


INTRODUCTION
The medical ultrasound is a diagnostic tool to provide images of organs and structures non invasively. Compared to other modalities such as CT scans and MRI, the ultrasound has the lower price, easy to move (portable), small and easily manipulated, and non -ionizing [1], [2]. Besides being used to monitor pregnancy, the ultrasonography is widely used for examination of diseases such as liver cancer [3], cardiovascular disease [4], breast cancer [5], [6], and musculoskeletal disorders [7]. Nevertheless, there are many drawbacks of this modality such as the artefacts on ultrasound images that named speckle noise. The presence of the noise will lead to the difficulty of further analysis of ultrasound images.
According to [8], a segmentation technique to obtain the fetal biometric measurement automatically on an ultrasound image could improve the efficiency of the clinical workflow. The abdominal circumference (AC) is fetal biometric size that used to confirm the gestational age, to estimate the fetal weight, and to identify the growth patterns and the abnormalities of the fetus. The existing studies of the AC segmentation technique included the edge detector approach, clustering technique, Hough transform,and gradient vector field snake (GVF) method [9]- [11]. In this work, the active geometric shape model (AGSM) approach that proposed by [1] was applied to segment the AC automatically. The AGSM could identify various geometric shapes of lines, circles, ellipse and cubic splines. The advantage of the AGSM method was robust to noise, broken lines and the non-uniform contrast [12]. The fetal abdomen image has no clear boundaries and has inconsistencies in the internal structure [13]. Moreover, an ultrasound image often contains speckle noise. So it takes a preprocessing technique called the speckle noise reduction method that aims to elimin ate noise while maintaining the feature of the boundary.
The speckle appears as a grain pattern that varies depending on the type of biological tissue. According to Joel & Sivakumar [14], there were two groups of the speckle noise reduction methods, i.e. the spatial filtering method and the multiscale or multiresolution method. Loizou & Pattichis [15] divided the spatial filtering techniques into three groups, i.e. the linear filtering approach such as the Lee filter [16] and the Frost filter [17], the nonlinear filtering technique such as the bilateral filter [18] and the shock filter [3], and the diffusion filtering method. The diffusion filtering techniques, such as the conventional anisotropic diffusion method [19] and the speckle reduction anisotropic diffusion (SRAD) method [20], adopted the partial differential equations and set up an adaptive filtering constraint in the denoising process. The SRAD approachwas the development of the anisotropic diffusion method. The last group of the speckle reduction methods was the multiresolution based method that one of which implements a wavelet transform. The wavelet based technique proposed by Babu et al [21] was used to remove the noise speckle from the SAR image. Whereas, Gopinathan et al [22] suggested the wavelet based denoising scheme to eliminate Gaussian noise. Hamiane & Saeed [23] applied the discrete wavelet transform to reduce noise in MRI brain images.
From all of the denoising methods, the SRAD is the most recent and widely developed technique. The SRAD is not only able to maintain the edge of the image but also able to increase the detail image. However, the SRAD is not able to eliminate the noise in a high-intensity area. Some hybrid techniques were proposed to overcome the denoising problem in the h igh-intensity area of an image. One of the hybrid denoising techniques was the multiresolution bilateral filter which combining the wavelet thresholding approach and bilateral filtering method [24], [25]. Zhang et al [25] applied the fast bilateral filtering to decrease the computational time of the technique proposed by Zhang & Gunturk [24].
In this paper, we propose a hybrid denoising framework that combining the multiresolution bilateral filtering method and the anisotropic diffusion method. To reduce speckle noise in the high-intensity area, we apply wavelet decomposition method that divides the image signal into different resolutions, i.e. an approximation subband and three detail subbands which represent the high -intensity area. We implement the anisotropic diffusion method in the detail resolution images to reduce noise better. This work does not only propose a speckle noise reduction method but also compare it to the other speckle noise reduction methods. This research also shows the influence of the proposed noise reduction technique in increasing the accuracy of the AC segmentation in the fetal ultrasound images.
The paper is arranged into four sections. In Section 1, we explain the problem that will be resolved and the literature review about the speckle noise reduction methods. Section 2 describes the research method that includes the existing denoising methods, the proposed denoising scheme, the image dataset, the AC segmentation method and the performance measurement. In Section 3, we analyze the experiment result, i.e. the noise reduction test and the segmentation test. Finally, Section 4 includes the conclusion of this work.

RESEARCH METHOD
We conducted two experiments to evaluate the proposed speckle noise reduction method as follows: (1) noise reduction experiment and (2) abdominal circumference segmentation experiment. The objective of the first test is to evaluate the capability of the proposed method in reducing the speckle noise and maintaining the edge. The second trial aims to show the influence of noise reduction techniques in the accuracy improvement of the abdominal circumference segmentation on the ultrasound image. In this section, we present a research method that it includes the existing despeckling techniques, the data, the proposed method, the fetal abdomen segmentation, and the performance evaluation technique .

Existing denoising methods
Firstly, we introduce the existing despeckling methods that are compared to the proposed approach as follows : a. Lee filter [16] Lee filter aimsto eliminatethe noise in a radar image while maintaining edges or features. This filter isbased on the linear filter model and the least mean square error approximation. Since, I x,y is the pixelintensityat the coordinates (x, y) in a MxN noisy image and ℵ x,y I x,y are pixels in a n x n neighborhood window, and m s is the mean of the intensities in the window ℵ, the Lee filter is expressed in the equation as follows : where d p,q is the distance from pixelpat coordinate(i p ,j p ) to pixel qat position (i q ,j q ), with p,qℵ x,y , ℵ x,y is the nxn neighbor window, that is expressed in the following equation: The variable K is the damping factor chosen such that the homogenous area KC s 2 close to zero.
c. Bilateral filtering [18] The bilateral filtering is a nonlinear filtering technique that can eliminate noise while preserving the edge of the image. Moreover, this filter implements the average of image intensities in a neighborhood window in a similar way to Gaussian filtering. To maintain the edge during the smoothing process, the bilateral method calculates the difference in a pixel intensity with its neighboring pixels. The following equation is used to express the bilateral function.
where I(p) is an input image and ℵ(p) is a local window centered at the pixel p,while qℵis the other pixel in the window ℵ. The function G σ (.) is a Gaussian function with σ as its standard deviation. The value and are parameters controlling the amount of the filtering in the image domains, and C is a constant to normalize the function such that the sum of weights equal to 1.The control parameters for denoising applications are determined using a linear relationship between and , where is a noise standard deviation.
d. Shock filter [3] Shock filter is based on a morphological operator, dilation or erosion, applied locally. The dilation or the erosion application dependson whether the pixel belongs . The influence zones, the ma ximum or minimum zone, aredetermined by the sign function. The sign function called the signum, s, in set {-1, 0, +1} is based on the Laplace operator. The border line between two influence zones produces a sharp discontinuity called shock.Shock filters applythe partial differential equation that run repeatedly as n iterations. Leta noisy imagef(x,y) and t=1,2,...,n. Thefiltered image equation at the t-th iterationis expressed as follows : (6) where = (u x , u y ) t is the spatial gradient of the function uat vertical direction (u x ) and at horizontal direction (u y ), and √ . The function is modified using the coherenceenhancing shock filtersthat is proposed by [26], such that the Equation (6) becomes: where v = K σ * u, and w is the normalized dominant eigen vector of the structure tensor * +, and v xx , v xy , v yy are the second order derivatives of v.
e. Anisotropic diffusion [19] and SRAD [20] The SRAD method is an extension of the anisotropic diffusion approach. The a nisotropic diffusion method uses the partial differential equation (PDE) approach which is a nonlinear filtering method that is expressed by the following equation: (9) where n is the number of iterations, f(i, j) is an original image functionin the (i, j) coordinates and c (x, y, z, t) is the diffusion coefficient with the horizontal, vertical and diagonal directionat the tth iteration. The diffusion coefficient allows to encourage the homogeneous image areasdiffusion and prevent the diffusion at the edges. The following equation expres ses the diffusion function that is used to enhance the higher contrast areas than, the lower ones.
Conversely, the second diffusion function that affects thelarger areas is expressed as follows: where k is a gradient threshold which controls the conduction. In the SRAD method, the k-value is estimated by the intensity average of a chosen homogenous area.At each step, the PDE is updated by the following equation : with is a integral constant, .
f. Decomposition wavelet [27] According to Sendur and Selesnick [27], there arethree main steps for the denoising algorithms using wavelet transforms, as shown in Figure 1. Firstly, the noisy image is analysed with discrete wavelet transforms (DWT). In this step, an input image is decomposed into the approximation subband using the lowpass filter (L) and the detail subbands using the highpass filter (H). The detail subbands provide the detailed information in the vertical (LH), horizontal (HL) and diagonal (HH) directions of the original image. Furthermore, the wavelet coefficients in the detail subbands are modified using thresholding processes (T). The last step is the synthesis process using the inverse discrete wavelet transform (IDWT). The wavelet thresholding aims to eliminate the wavelet coefficients of the detail subbands that are lower than a particular threshold value. This process impacts on noise reduction without affecting the other essential features in the image. There are two thresholding methods including the soft thresholding and hard thresholding. In the soft thresholding procedure that is proposed by Donoho [28], the wavelet coefficients that smaller than the threshold value are set to zero, otherwise, the coefficients are subtracted from the threshold value. The following equation expressesthe soft thresholding.
The threshold value T is estimated by the Bi variate shrinkage [27] using the following formula: where w i is the diagonal wavelet coefficientsof the ith decomposition level.
g. Multiresolution bilateral filtering [24] The multiresolution bilateral filteringis a hybrid method that integrates the multiresolution wavelet decomposition and bilateral filtering method. The bilateral filtering is appliedto the approximation subband of wavelet decomposition. Meanwhile, the detail sub-bands of wavelet decomposition are computed by the wavelet thresholding. The threshold value is calculated by the estimation threshold functionas shown in Equation (14). In their paper, Zhang & Gunturk [24] identifiesthe linear relation between the characteristic of intensity domain and the noise standard deviation . The optimal value of the parameter is related to with scale factor sas expressed by the following equation: (15) Meanwhile, is defined as a ratio between the median absolute deviation (MAD) of the detail subband (HH) coefficients and a decimal number 0.6745.

Data
In this research, we used two kinds of imagesfor each experiment. The first type was a phantom image with the size of 200x200, the modified Shepp-Logan, that was generated by the phantom function in Matlab 2013 as shown in Figure 2(a). This phantom image was added with speckle noise for noise variance equal to 0.1, as shown in Figure 2(b). The other image was the real abdomen ultrasound image that was obtained using a portable scanner, i.e. the SonoStar UBox-10 Ultrasound B Scanner with a transabdominal convex probe 3.5 MHz. The example of the real fetal abdomen image is presented in Figure 3.

Proposed speckle noise reduction method
The proposed method for ultrasound image denoising described in Figure 4 includes five processes: the discrete wavelet transform (DWT), bilateral filtering, the thresholding process followed by anisotropic diffusion, and the inverse wavelet transform (IDWT). Firstly, the noisy image is decomposed into the nth level by the discrete wavelet transform (DWT) that generates a subband which is an approximation of the original image and the three sub-bands that provides detailed information on the vertical, horizontal and diagonal of the original image. Mother wavelet function applied in this method is a Daubechies "db8" function [29]. The approximation subband is processed by bilateral filtering in the Equation (5). The detail subbands are computed using the soft thresholding technique as in Equation (13) and followed by the anisotropic diffusion method with the second diffusion function as in Equation (11). To estimate thethreshold value, we use the formula in Equation (14). Unlike other existing anisotropic diffusion in [19] and [20], we apply the derivative filters for eight directions at the derivative approximationstage of the anisotropic diffusion process. The anisotropic diffusion process can increase noise removal in detail subbands. After all the process, the inverse wavelet transform is applied to get an approximation image. All of these steps arerepeated using the approximation image as the input until it reaches the first level decomposition wavelet. Figure 4 shows an example of denoising image for level n=1. To achieve the best result, we determine the parameters required in this method as shown in Table 1.

Fetal abdominal circumference segmentation method
The fetal abdomen segmentation in an ultrasound image is used to get the abdominal circumference (AC) with standard measurement as presented in Figure 5 where the left figure shows the landmark of the abdomen image in the right figure. The landmark contains (1) umbilical vein (2) middle portal vein (3) stomach (4) spine (5) rib. The abdominal circumference is measured around the outside of the skin line as shown by white circle in Figure 5 [30].
We apply the active geometric shape model (AGSM) method that is proposed by Wang & Boyer [12] to detect and to measure the fetal abdominal circumference. The AGSM method implements a geometric shape approach represented by parametric equations to fit an object in an image. This technique adjusts the shape parameters according to the integral of a force field along the shape contour. The output of the AGSM function is the radius of a generated circle that applied to calculate the circumference. Figure 6 shows the block diagram of the abdominal circumference segmentation process. The input images are the ultrasound image from denoising process that cropped automatically. At the first step, the Canny edge detector is applied to get the boundaries in the binary image and afterwards the result is processed by the Gaussian smoothing filter. Furthermore, the gradient vector field (GVF) is derivedfrom the blurred image in the vertical direction and horizontal direction. Based on the two vector fields and initial circle parameters, a circle equation is updated.

Quantitative quality measurements
In this research, we used four quality performance in the first experiment, i.e., the peak signal to noise ratio (PSNR) [3], the root mean square error (RMSE) [3], the Pratt"s Figure Of Merit (FOM) and the Mean Structural Similarity Index Measure (MSSIM). The MSSIM that proposed by [31] is a measurement used to measure the image restoration results quality that corresponds to the human visual perception by combining three common factors that are the loss of correlation, luminance and contrast distortion. The Pratt"s figure of merit (FOM) that submitted by [20] is a measurement to compare edge preservation performance of different speckle reduction schemes. Results of FOM depend on the edge detection method used. In this research, we implemented the Canny Edge Detector.
In the segmentation experiment, we measured the performance using the accuracyand the Dice similarityas used in [13]. Suppose denotes the result of segmentation for method M and is ground truth obtained from manual segmentation. The accuracyis the average of two values , i.e., True Positive (TP) and True Negative (TN). TP is formulated as the ratio of the intersection between and to the region covered by . TN is the ratio of the region not covered by and to the areas outside the region . True-Positive Ratio (TPR) also known as sensitivity and True-Negative Ratio (TNR) also known as specificity [6]. Dice similarity (D) is the primary similarity metric used to measure overlapping between two regions spatially, defined as follows.

Reduction noise performance
In the first experiment, we compared our method with existing methods using the parameters that contained in Table 2 with the PSNR, RMSE, MSSIM, and Pratt's FOM quantitative measurements. To obtainthe Pratt"s FOM, we applied the Canny edge detector with the standard deviation of the Gaussian filter σ=√ to the original phantom image and noisy phantom image. Table 3 summarizes the comparison of the quality measurement results from various methods. The bilateral filtering method has the best value of the mean square error and the peak signal to ratio values. Unfortunately, from the Pratt's FOMvalue, the bilateral filterability to maintain the edge is lower than other methods. Nevertheless, usage of bilateral filtering will be able to maintain the visu al structure that corresponds to human visual perception, as seen from the MSSIM. The thresholding approach applied to the multiresolution wavelet based method can keep the average error small in the image restoration process, as seen from the RMSE value. Therefore, the integration of the bilateral filter, the anisotropic diffusion and the thresholding wavelet in our proposed method offer the best performance quantitatively. T he MSSIM index value of the proposed method gives a significant improvement. Consequently, the proposed method result has a great similarity structurallyto the reference image.Based on the Pratt's FOM, this method provides a better value. Therefore, it indicates that the method is able to preserve the edge. Frost filter result has the high Pratt"s FOM value, conversely, the PSNR and the RMSE value of the Frost filter methodare worse than the proposed method. Visually in Figure 7, the proposed method gives the best image that can reduce noise and remain preserved the edges. In Figure 7(c), Frost filter result looks more blurred, such that the features or the edges of the image are lost.

CONCLUSION
The proposed denoising method was a hybrid speckle noise reduction framework that combinedthe anisotropic diffusion method and multiresolution bilateral filtering method. In this proposed scheme, the wavelet thresholding technique and anisotropic diffusion method were applied to the detail subbands of the wavelet decomposition. Meanwhile, the approximation subband was processed by the bilateral filtering method. To find the performance of the proposed approach, we compared it with the other denoising techniques, such as the Lee filter, Frost filter, shock filter, bilateral filter, wavelet thresholding, and multiresolution bilateral filtering method. The comparison was made regarding noise reduction and the effect of noise reduction results to improve the abdominal circumference segmentation performance in fetal ultrasound images. From Pratt"s FOM measurement, this method obtained the best performance to remove speckle noise without reducing the edge. Moreover, this approach also had the best quality of human visual perception. For the segmentation task, the proposed denoising method could improve the accuracy of the fetal biometric abdominal circumference segmentation. The segmentation performance variance of the intended image is quite small, so it means that the proposed method can enhance the various types and quality of ultrasound images. For further development, an attempt to integrate other wavelet functions, such as the complex wavelet transform, rational dilation wavelet transform, or stationary wavelet transform.