Multiresolution generalized N dimension PCA for ultrasound image denoising

Background Ultrasound images are usually affected by speckle noise, which is a type of random multiplicative noise. Thus, reducing speckle and improving image visual quality are vital to obtaining better diagnosis. Method In this paper, a novel noise reduction method for medical ultrasound images, called multiresolution generalized N dimension PCA (MR-GND-PCA), is presented. In this method, the Gaussian pyramid and multiscale image stacks on each level are built first. GND-PCA as a multilinear subspace learning method is used for denoising. Each level is combined to achieve the final denoised image based on Laplacian pyramids. Results The proposed method is tested with synthetically speckled and real ultrasound images, and quality evaluation metrics, including MSE, SNR and PSNR, are used to evaluate its performance. Conclusion Experimental results show that the proposed method achieved the lowest noise interference and improved image quality by reducing noise and preserving the structure. Our method is also robust for the image with a much higher level of speckle noise. For clinical images, the results show that MR-GND-PCA can reduce speckle and preserve resolvable details.


Introduction
An accurate anatomical display of organs plays a very important role in current-day medical diagnosis. Among a diverse set of medical imaging modalities, the use of ultrasound has increased rapidly because of its non-invasive nature, hardware portability, inexpensiveness, and real-time imaging. However, ultrasound images have much more noise than other medical imaging modalities, such as computer tomography (CT) and magnetic resonance imaging (MRI). The primary source of ultrasound imaging noise is considered as speckle noise, which deteriorates image quality, degrades fine details and edge definition, and limits contrast resolution [1]. Consequently, reducing speckle noise while preserving anatomic information is necessary to reliably delineate the regions of interest and increase the diagnostic potential of ultrasounds.
Extensive work using various techniques has been conducted to reduce speckle noise. These techniques can be classified into two categories according to the moment when speckle reduction is produced, namely, compounding approach and post-processing approach [2]. Compounding approach modifies the data acquisition procedure to are invariant to changes in the spatial scale. Additionally, the larger scale structure is lost when the image is decomposed into small patches. This case can be avoided by employing the multiresolution/multiscale approach. Instead of traditional PCA, a multilinear subspace learning method called GND-PCA is also utilized to preserve useful information and noise suppression. The proposed method has three stages. First, a Gaussian pyramid for the input noisy image is built with a multiresolution scheme. The noise in the ultrasound image is gradually smoothed and computational time is reduced when computing in higher pyramid levels. The stacks are then built with multiscale images produced from each pyramid level. Second, dense patches (two-way array) with the same size are extracted, and similar patches are grouped together in three-way arrays in different stacks. GND-PCA is directly used to filter the three-dimensional array with its original structure in order to overcome the drawbacks of conventional PCA, as well as to consider the interrelationship among different modes of each patch to maintain useful information and remove noise. Denoised patches are then placed back into their original positions. Third, after obtaining the representative images in each stack, a combination motivated by Laplacian pyramids is obtained. And then, the denoising image is yielded from the combined image with exponential transformation and sharpen procedure.
The organization of the paper is as follows. Section Related work provides a brief introduction of related work. The proposed MR-GND-PCA denoising method is outlined in Section Methodology. The experimental design is presented and results are discussed in Section Experiments and results. Finally, Section Discussion and conclusions concludes the work.

Related work
Generalized N dimension PCA (GND-PCA) is the main related work of our proposed method. For more details please refer to [18]. The GND-PCA is briefly presented as follows.
Traditional PCA [19] is one of the fundamental linear subspace learning techniques for data analysis. The main purpose of PCA is to reduce data dimension and retain information that characterizes the variation of data as much as possible. Thus, PCA can be widely used in noise reduction, i.e., to keep useful information and eliminate high-frequency noise. However, the data must be unfolded into one-dimensional vectors to fit with the PCA processing. For a two-dimensional image or N-way array, this unfolding process destroys the structure of the input data. PCA also suffers from a generalization problem because the calculated bases with PCA cannot represent data very well if the number of samples is much smaller than the dimension of the unfolded vector.
Xu et al. [18] proposed a multilinear subspace learning method called generalized N dimension PCA (GND-PCA) to overcome these insufficiencies. GND-PCA is proposed by extending the conventional linear PCA based on multilinear algebra. In multilinear algebra, higher-dimensional data are treated as higher-order tensors that are N-way array [20]. The order of a tensor is equal to the number of ways, which are known as modes. GND-PCA directly analyzes the N-way array on the sub-mode spaces. This step not only retains the structure of the data, but also overcomes the generalization problem that occurs when the dimension of the unfolded data is much larger than the number of samples.
The bases are calculated on each mode subspace in order to provide a compact representation of the tensors. An iteration algorithm is utilized to obtain the optimal base matrices for each mode. Next, the core tensor containing all principal components is attained with the optimal bases. In this paper, we denote scalars by lower-case letters (x, y, …), vectors (one-way array) by boldface letters (x, y, …), matrices (two-way array) by boldface capital letters (X, Y, …), and third-order tensor (three-way array) by calligraphic capital letters (X ; Y; …).
GND-PCA was originally proposed for the construction of statistical appearance models. According to the intrinsic algorithm that lower-rank tensors are found to approximate the original tensors with more compact and meaningful form, we can utilize GND-PCA to retain useful information while suppressing redundancy, that is, image denoising. The next section will introduce our proposed noise reduction method in detail.

Methodology
Motivated in part by image pyramids and GND-PCA, we propose a novel ultrasound image denoising method, called multiresolution Generalized N Dimension PCA (MR-GND-PCA). Figure 1 shows the proposed MR-GND-PCA algorithm with its three stages.
(1) In the multiresolution decomposition stage, logarithmical transformation is first used to convert the multiplicative speckle noise to an additive noise. Then, a Gaussian pyramid for the noisy image is built using the multiresolution scheme. Multiscale image stacks on each level are built for each Gaussian pyramid. (2) In the noise reduction stage, similar patches (two-way array) are extracted in each stack and grouped into several categories (three-way array). GND-PCA is directly utilized for three-way arrays to suppress the noise. The denoised patches are put back in the original places to obtain clean multiscale image stacks. (3) In the multiresolution combination stage, the representative image on each level can be accessed by image averaging. The combination of representative images followed by exponential transformation and sharpen procedure will yield the noise reduction ultrasound image.
The above stages will be formulated in detail in the following subsection.

Multiresulution decomposition
A reliable model of speckle noise is presented in [21,22] g or x; y ð Þ ¼ S or x; y ð Þη or mu x; y ð Þþη or ad x; y ð Þ ð1Þ where S or (x, y) is a noise-free original gray-valued image, g or (x, y) is an observed noisy gray-valued image, η or mu x; y ð Þ is multiplicative noise (coherent interfering), η or ad x; y ð Þ is additive noise (such as sensor noise), and x and y are variables of spatial locations, (x, y) ∝ R 2 . We could approximate the Eq. 1 as Eq. 2, since the effect of additive noise is much smaller than that of multiplicative noise.
g or x; y ð Þ≈S or x; y ð Þη or mu x; y ð Þ ð2Þ Logarithmic transformation can be executed on both sides of Eq. 2 to separate the noise from the original image.
We simplify Eq. 3 as Eq. 4 without modifying the definition: Therefore, multiplicative speckle noise is converted to additive noise with logarithmical transformation. This step is more reasonable for the proposed MR-GND-PCA as applied to ultrasound image denoising. The Gaussian pyramid for the noisy image in the logarithmic space and the multiscale image stacks on each pyramid level are built successively because only the possible scale-space kernel is a Gaussian function [23]. The scale space of the image g is produced from the convolution of variable-scale Gaussian filters f(σ, (x, y)) with the input image: where "*" is the convolution operation in x and y, and σ denotes the scale parameter, and In our paper, an associated Gaussian pyramid G a,0 , 0 ≤ a ≤ A is composed of a set of different resolutions of g. Here, A + 1 is the number of total levels and G 0,0 = g. G a,0 are created by applying the Gaussian filter and successive downsampling. Both operations are combined in the "Reduce" operation, and shown on the first line in Figure 2, where f 1 denotes a linear low-pass Gaussian filter, and (2 ↓) denotes downsampling that discards every data element with an odd index and is applied to each dimension of the image. Multiscale image stacks on each pyramid level are subsequently constructed. Assume G a,0 as the basic image of the a th stack at each level of the Gaussian pyramid. G a,b , 0 ≤ b ≤ B are a series of different corresponding scales of G a,0 ; B + 1 is the number of scales on each level. The multiscale image stacks G a,b , 0 ≤ b ≤ B are created by applying the low-pass Gaussian filter f 2 . Thus, we have where (Num) denotes the number of f 2 convolutions with G a,0 . Both Gaussian filters f 1 and f 2 are defined in 5 × 5 sliding windows.

Noise reduction
GND-PCA as a multilinear subspace method is utilized for denoising the images of each stack on the pyramid level to reduce noise for the entire image. Take the first stack, G 0,b , 0 ≤ b ≤ B, for example. Q overlapping dense patches, X 0 ð Þ q ∈R I 1 ÂI 2 ; q ¼ 1; 2; …; Q, are extracted from all images in this stack. Here, I 1 × I 2 is the size of the patch, and the superiors "(0)" denotes that all the extracted patches are in the first stack. Organizing the patches into sensible groupings is a fundamental process for GND-PCA denoising. Cluster analysis, which consists of hard and soft clustering, implements groupings according to measured intrinsic characteristics or similarities. Loosely speaking, hard clustering assigns the data point to one and only one cluster. In contrast, soft clustering assigns the data point to all clusters with different weights whose sum equals one. Clustering is unsupervised learning due to the absence of category information. Despite the fact that thousands of clustering algorithms have been published, we adopt the most popular and simple clustering algorithm, K-means [24], for our research. As the discussion of diversified clustering algorithms is beyond our scope, please refer to [25] for more details on data clustering.
The advantage of containing the original measured quantities and easy implementation allows for pixel intensity values instead of patch features to be employed for clustering. The overlapping dense patches on the first stack are unfolded as the set of ; Q that will be clustered into K clusters. Nevertheless, massive noise existing in the patches decreases the reliability of the characteristic structure and eliminates the performance of the clustering algorithm. In consideration of computational effectiveness, (I 1 × I 2 ) -dimensional points are projected into the lower dimensional space by PCA to characterize the intrinsic variation of the data. After projection, ( denotes the mean of cluster c k . A partition of the patches can be found by minimizing the squared error between the grouping means and the points: After clustering in one multiscale image stack, the corresponding patches that exhibit high correlation are grouped into K three-way arrays X 0 ð Þ i ∈R I 1 ÂI 2 ÂI 3 ; i ¼ 1; 2; …; K , with I 3 denoting the number of similar images in the group. Three-way arrays are processed precisely by GND-PCA to suppress noise without destroying the internal connection of patch similarities.
A three-way array X 0 ð Þ i ∈R I 1 ÂI 2 ÂI 3 ; i ¼ 1; 2; …; K can be unfolded into three matrices along each mode space and represented by X , where, n = 1, 2, 3. Generally, the number of D i(n) are decided by the experience and application. A simple determination procedure proposed in [26] can also be used. In each mode, this method decides the dimension by the variation ratio to be kept by eigenvalues. The order of the eigenvectors, u Dimension reduction is completed in each mode for denoising, and thus, the calculated U 0 ð Þ i n ð Þ are not the optimal bases for the entire three-way array. The optimal bases for the three-way array and denoised patches can be obtained by utilizing the above U 0 ð Þ i n ð Þ as the initialization. Then, an iteration algorithm, which calculates one matrix if the others are fixed, is used to determine all optimal matrices to satisfy the following energy function: Þ ÂI 3 are core tensors. Note that the dimension of the third mode is retained without reduction. The reconstructed three-way array shown in Eq. 11 is the denoised grouping. Figure 4 illustrates the dimension reduction of GND-PCA with a three-way array.
After all the three-way arrays on the same level are denoised with GND-PCA, we place them back into their original positions. With the overlapping patches extracted, more than one patch estimation in the denoised three-way arrays can be located at exactly the same coordinate. An accumulation is performed by averaging the overlapping pixel positions. Take the example on G 0,0 , which denotes the first level image with no Gaussian filter. The denoised estimationĜ 0;0 of G 0,0 is computed by an average of the patch estimatesŶ 0 ð Þ p ∈R I 1 ÂI 2 ; p ¼ 1; 2; …; P in the first stack. Here,Ŷ 0 ð Þ p ∈Ŷ 0 ð Þ q ; P≤Q. We havê Figure 4 Dimension reduction of GND-PCA with a three-way array. where E has the same size asŶ 0 ð Þ p and all elements of E are one, and x R is the coordinate iterated in patches. The denoised estimations for all stacks are denoted bŷ G a;b ; 0≤a≤A; 0≤b≤B.

Multiresolution combination
All images on the same stack are added up to become one representative image, as shown in Eq. 13.
The denoised image is obtained by combining the representative images G a , a = 0, 1, 2, …, A at each level. The corresponding pyramids with levels L a contain the differences of high-frequency components between G a and G a + 1 . Differences G a − G a + 1 , a = 0, 1, …, A − 1 cannot be calculated directly because G a + 1 contains fewer samples than G a . For this reason, the number of samples in G a + 1 is increased to match G a by interpolating the missing samples and the interpolation filter. Both operations are combined in the "Expand" operation, where (2 ↑) denotes upsampling by inserting a zero between adjacent data elements and f 3 denotes the linear interpolation filter. The output imageĝ can be reconstructed from L a by, and then, In the reconstructed process, the low-frequency information of G a is discarded and replaced by G a + 1 . Finally, the de-speckled image in log-domain is converted into spatial domain through exponential transformation. The complete MR-GND-PCA method is shown as a diagrammatic sketch in Figure 2.

Experiments and results
The proposed MR-GND-PCA was used as an application to denoise an ultrasound image corrupted by speckle noise. The quantitative performance of the proposed method was investigated by using three different kinds of images, including a synthetic image with artificial speckle noise, phantom data, and a common carotid artery scanned by commercial ultrasonic scanners.

Synthetic image investigation
The efficiency of the proposed MR-GND-PCA was tested by using a generated noisy image with artificial speckle noise. Equation 2 shows that we degraded the original test image by multiplying it with unit mean random fields to obtain the test speckle image. The spatially correlated speckle noise was produced by low-pass filtering a complex Gaussian random field and taking the magnitude of the filtered output [27]. Low-pass filtering was performed by averaging the complex values of the size three sliding window. This short-term correlation allowed the noise correlation to taper gradually to zero, as well as the correlation length of the speckle to be controlled by appropriately setting the size of the kernel used to introduce correlation to the underlying Gaussian noise. The size of the synthetic image is 141 × 141. Two distinct levels of speckle noise were generated by changing the variance of the underlying complex Gaussian random field. The purely synthetic image shown in Figure 5 is invoked as a reference noise-free image and is available at [34]. The image is comprised of regions with uniform intensity, sharp edges, and strong scatters. The performance of the proposed MR-GND-PCA method was compared with traditional despeckling methods (Median filter, Fourier ideal filtering, Fourier Butterworth filter, Wavelet filter) as well as state-of-the-art algorithms (Bilateral filter [28], KSVD [12], Non-local means [29], BM3DPCA [13], PGPCA [14], LPGPCA [15]). All despeckling methods were implemented on the images with high, middle and low levels of speckle noise.
We used four well-known metrics as quantitative performance measures to assess quality image, including mean-square error (MSE), signal-to-noise ratio (SNR), peak signal-to-noise ratio (PSNR), and the structural similarity (SSIM). MSE evaluates the noise variance and is given by whereĝ x; y ð Þ is the denoised image and S(x, y) is the purely synthetic image. SNR is defined as where P s is the variance of the noise-free reference image as PSNR is defined by where N max is the maximum fluctuation in the input image.  And SSIM is given by where μĝ x;y ð Þ , μ S(x,y) , σĝ x;y ð Þ , σ S(x,y) , and σĝ x;y ð ÞS x;y ð Þ are the mean ofĝ x; y ð Þ, the mean of S (x, y), the variance ofĝ x; y ð Þ, the variance of S(x, y), and the covariance ofĝ x; y ð Þ and S (x, y), respectively; C 1 and C 2 are small constants.   , which are also considered as the best parameters in the following experiments, though SSIM is a little bit higher with σ = 1 than that with Figure 9.

Comparison of MR-PCA and MR-GND-PCA
In order to illustrate the effectiveness of GND-PCA for image denoising, we substitute GND-PCA for PCA in the proposed MR-GND-PCA framework, which is called MR-PCA. The restoration results for low-, middle-, and high-level noise are shown in Figure 10.   MR-GND-PCA gives better overall performance in terms of both noise suppression and detail preservation. However, MR-PCA produces artifacts around the edges of each objects, which is a potential influence of the accurate diagnosis.

Comparison of multiple denoising methods
The values of MSE, SNR, PSNR and SSIM obtained for the corresponding test image with all the methods are summed up in Tables 1, 2 and 3.  Tables 1, 2  the denoising methods, such as upper right object, the proposed MR-GND-PCA is the most reliable process algorithm because of suppressing the majority noise keeping the details in the source image. The result confirms that the multiresolution created for GND-PCA performance is an appropriate combination for ultrasound images with massive noise compare to the state-of-art methods.

Phantom data validation
Phantom data called CIRS MODEL 057A are available at http://www.cirsinc.com/. The B-scan image of these phantom data was obtained from a commercial ultrasonic scanner with a 3.

Clinical data verification
Clinical images of common carotid artery (CCA) are available for free on the Signal Processing Laboratory website. The B-scans were acquired with a Sonix OP ultrasound scanner with a different set-up of depth, gain, time gain compensation (TGC) curve, and different linear array transducers. Figure 15 shows one of the CCA images with size of 330 × 528, in which three highlighted areas marked by the yellow arrows can be fully identified and preserved. Figure 16 illustrates the results processed before and after using MR-GND-PCA. The highlighted areas in Figures 16(c) and (d) show that the denoised image in Figure 16(b) displays significant improvement and recovers discontinuities well.

Discussion and conclusions
In this paper, we introduced MR-GND-PCA for noise suppression in ultrasound images. Our approach consists of three stages. First, after the ultrasound image is converted by logarithmical transformation, the Gaussian pyramid with multiscale image stacks on each level is built. Second, similar patches are grouped to be directly denoised by GND-PCA in the same pyramid level and then placed on the original location to form clean pyramid levels. Third, the representative images obtained from each level are combined to obtain the noise reduction ultrasound image. The first and third stages create a multiresolution framework that utilizes all frequency bands. Downscaling has a denoising effect that contributes to better denoising. The second stage uses GND-PCA to maintain the structure of the data and to consider the correlation of all modes of the proposed data.  Experiments were carried out on synthetic images with three distinctive level speckles and two clinical images. Quantitative measures were used to compare 11 denoising filters for the synthetic images. Experimental results show that our proposed method achieved the lowest noise interference. Our method is also robust for the image with a much higher level of speckle noise. For clinical images, the results show that MR-GND-PCA can reduce speckle and preserve resolvable details. Further work on speeding up the method and enhancing the edge of the objects will be pursued in the future.