 Research
 Open Access
 Published:
Multiresolution generalized N dimension PCA for ultrasound image denoising
BioMedical Engineering OnLinevolume 13, Article number: 112 (2014)
Abstract
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 (MRGNDPCA), is presented. In this method, the Gaussian pyramid and multiscale image stacks on each level are built first. GNDPCA 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 MRGNDPCA can reduce speckle and preserve resolvable details.
Introduction
An accurate anatomical display of organs plays a very important role in currentday medical diagnosis. Among a diverse set of medical imaging modalities, the use of ultrasound has increased rapidly because of its noninvasive nature, hardware portability, inexpensiveness, and realtime 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 postprocessing approach [2]. Compounding approach modifies the data acquisition procedure to produce several images of the same region and then combine them to form a single image. This approach is beyond the scope of this paper. Postprocessing approach processes ultrasound images after being generated. Most studies have focused on this approach, in which filtering methods are developed to be applied directly to the original image. Based on the assumption that speckle is essentially a multiplicative noise [3], filtering methods composed of many fixed and adaptive filters have been proposed, such as L2mean filter [4], adaptive filter reduction [5], adaptive weighted median filter [6], nonlinear diffusion [7], Map estimation [8], and so on.
Although techniques for processing multiplicative noise exist, they are not as well developed as the techniques for additive noise. Generally, the multiplicative nature of speckle noise formation can be converted into an additive noise through the application of logarithmic transformation to a speckled image. Subsequently, various types of stateoftheart denoising algorithms can be utilized for ultrasound despeckling. Stateoftheart algorithms are mainly divided into three categories: (1) algorithms applied in the frequency domain, (2) algorithms applied in the spatial domain, and (3) algorithms applied in both spatial and frequency domains. Waveletbased filters implement denoising in the frequency domain, such as Bayesian least squares (BLS) Gaussian scald mixture (GSM) wavelet denoising method [9]. This method transforms the image into the wavelet domain, assumes the GSM model on neighborhoods, and denoises using BLS estimation. However, a rich amount of different local structures are contained in medical images, and these structures cannot be well represented by only one fixed wavelet basis. Thus, numerous visual artifacts can be introduced through waveletbased methods. Using the temporal perspective, spatial denoising methods are proposed to discover redundancies in a single image, such as nonlocal meanbased methods [10, 11]. These methods extract the commonalities and use a weighted average based on similarity. Recently, denoising methods implemented in both spatial and frequency domains, such as KSVD method [12], nonlocal collaborative filtering [13], and principle component analysis (PCA)based methods [14, 15], have attracted attention. The core idea of these denoising methods is the use of sparser representation instead of a spatial weighted average [16, 17].
Recent studies have shown that PCAbased methods can obtain impressive achievements. Bags of patches extracted from the noisy image are processed by PCA to reduce the redundancy between data. PCA is a linear subspace learning method that has been applied widely for dimensionality reduction by searching for the maximum variance directions. However, conventional PCA was originally proposed to process onedimensional vectors, which require all input data to first be unfolded into a vector in order to fit the nature of the PCA. For similar patches extracted to reduce redundancy, vectorization destroys the structure of the input patches, resulting in the overfitting problem because the dimension of the vectorized data may be larger than that of the sample number.
This paper proposes a novel multiresolution generalized N dimension PCA for ultrasound image denoising called the MRGNDPCA. The issue of noise for ultrasound images can be reduced with higher frequency imaging, but this would limit the depth of ultrasound penetration. Thus, the noise is decided by the different penetrations. An intuitive justification is that ultrasound images are equally likely to be viewed from different distances. Our proposed method assumes that the statistics of ultrasound images 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 GNDPCA 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 (twoway array) with the same size are extracted, and similar patches are grouped together in threeway arrays in different stacks. GNDPCA is directly used to filter the threedimensional 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 MRGNDPCA 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 (GNDPCA) is the main related work of our proposed method. For more details please refer to [18]. The GNDPCA 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 highfrequency noise. However, the data must be unfolded into onedimensional vectors to fit with the PCA processing. For a twodimensional image or Nway 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 (GNDPCA) to overcome these insufficiencies. GNDPCA is proposed by extending the conventional linear PCA based on multilinear algebra. In multilinear algebra, higherdimensional data are treated as higherorder tensors that are Nway array [20]. The order of a tensor is equal to the number of ways, which are known as modes. GNDPCA directly analyzes the Nway array on the submode 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 lowercase letters (x, y, …), vectors (oneway array) by boldface letters (x, y, …), matrices (twoway array) by boldface capital letters (X, Y, …), and thirdorder tensor (threeway array) by calligraphic capital letters ($\mathcal{X},\mathcal{Y},\phantom{\rule{0.5em}{0ex}}\dots $).
GNDPCA was originally proposed for the construction of statistical appearance models. According to the intrinsic algorithm that lowerrank tensors are found to approximate the original tensors with more compact and meaningful form, we can utilize GNDPCA 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 GNDPCA, we propose a novel ultrasound image denoising method, called multiresolution Generalized N Dimension PCA (MRGNDPCA). Figure 1 shows the proposed MRGNDPCA 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 (twoway array) are extracted in each stack and grouped into several categories (threeway array). GNDPCA is directly utilized for threeway 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]
where S^{or}(x, y) is a noisefree original grayvalued image, g^{or}(x, y) is an observed noisy grayvalued image, ${\eta}_{\mathit{mu}}^{\mathit{or}}\left(x,y\right)$ is multiplicative noise (coherent interfering), ${\eta}_{\mathit{ad}}^{\mathit{or}}\left(x,y\right)$ 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.
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 MRGNDPCA 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 scalespace kernel is a Gaussian function [23]. The scale space of the image g is produced from the convolution of variablescale 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 lowpass 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 lowpass 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
GNDPCA 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, ${\mathbf{X}}_{q}^{\left(0\right)}\in {R}^{{I}_{1}\times {I}_{2}},q=1,2,\dots ,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 GNDPCA 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, Kmeans [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 (I _{1} × I _{2}) dimensional points $\left[{\widehat{\mathbf{x}}}_{1}^{\left(0\right)},{\widehat{\mathbf{x}}}_{2}^{\left(0\right)},\dots ,{\widehat{\mathbf{x}}}_{Q}^{\left(0\right)}\right],q=1,2,\dots ,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, (I _{1} × I _{2}) dimensional points $\left[{\widehat{\mathbf{x}}}_{1}^{\left(0\right)},{\widehat{\mathbf{x}}}_{2}^{\left(0\right)},\dots ,{\widehat{\mathbf{x}}}_{Q}^{\left(0\right)}\right]$ is converted into J dimensional points $\left[{\mathbf{x}}_{1}^{\left(0\right)},{\mathbf{x}}_{2}^{\left(0\right)},\dots ,{\mathbf{x}}_{Q}^{\left(0\right)}\right]$.
Assume ${\mathbf{C}}^{\left(0\right)}=\left[{\mathbf{c}}_{1}^{\left(0\right)},{\mathbf{c}}_{2}^{\left(0\right)},\dots ,{\mathbf{c}}_{k}^{\left(0\right)},\dots ,{\mathbf{c}}_{K}^{\left(0\right)}\right],k=1,2,\dots ,K$ denotes K clusters and ${\mu}_{k}^{\left(0\right)}$ 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 threeway arrays ${\mathcal{X}}_{i}^{\left(0\right)}\in {R}^{{I}_{1}\times {I}_{2}\times {I}_{3}},i=1,2,\dots ,K$, with I _{3} denoting the number of similar images in the group. Threeway arrays are processed precisely by GNDPCA to suppress noise without destroying the internal connection of patch similarities.
A threeway array ${\mathcal{X}}_{i}^{\left(0\right)}\in {R}^{{I}_{1}\times {I}_{2}\times {I}_{3}},i=1,2,\dots ,K$ can be unfolded into three matrices along each mode space and represented by ${\mathbf{X}}_{i\left(1\right)}^{\left(0\right)}\in {R}^{{I}_{1}\times \left({I}_{2}\cdot {I}_{3}\right)}$, ${\mathbf{X}}_{i\left(2\right)}^{\left(0\right)}\in {R}^{{I}_{2}\times \left({I}_{3}\cdot {I}_{1}\right)}$, and ${\mathbf{X}}_{i\left(3\right)}^{\left(0\right)}\in {R}^{{I}_{3}\times \left({I}_{1}\cdot {I}_{2}\right)},i=1,2,\dots ,K$. The unfolding procedures, called mode matricization of a thirdorder tensor, are visualized in Figure 3. For each unfolded matrix, the eigenvector matrix ${\mathbf{U}}_{i\left(n\right)}^{\left(0\right)}=\left[{\mathbf{u}}_{i\left(n\right)1}^{\left(0\right)},{\mathbf{u}}_{i\left(n\right)2}^{\left(0\right)},\cdots ,{\mathbf{u}}_{i\left(n\right){D}_{i\left(n\right)}}^{\left(0\right)}\right]$ associated with the first ${D}_{i\left(n\right)}^{\left(0\right)}$ largest eigenvalues (number of modesubspace bases) can be calculated from the covariance matrix ${\mathbf{C}}_{{\mathbf{X}}_{i\left(n\right)}^{\left(0\right)}}=E\left\{{\mathbf{X}}_{i\left(n\right)}^{\left(0\right)}{{\mathbf{X}}_{i\left(n\right)}^{\left(0\right)}}^{T}\right\}$, 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, ${\mathbf{u}}_{i\left(n\right)1}^{\left(0\right)},{\mathbf{u}}_{i\left(n\right)2}^{\left(0\right)},\cdots ,{\mathbf{u}}_{i\left(n\right){D}_{i\left(n\right)}}^{\left(0\right)}$, satisfies the corresponding eigenvalues ${\lambda}_{i\left(n\right)1}^{\left(0\right)}\ge {\lambda}_{i\left(n\right)2}^{\left(0\right)}\ge \dots \ge {\lambda}_{i\left(n\right){D}_{i\left(n\right)}}^{\left(0\right)}$. Dimension reduction is completed in each mode for denoising, and thus, the calculated ${\mathbf{U}}_{i\left(n\right)}^{\left(0\right)}$ are not the optimal bases for the entire threeway array. The optimal bases for the threeway array and denoised patches can be obtained by utilizing the above ${\mathbf{U}}_{i\left(n\right)}^{\left(0\right)}$ 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:
where ${\mathcal{Y}}_{i}^{\left(0\right)}\in {R}^{{D}_{i\left(1\right)}^{\left(0\right)}\times {D}_{i\left(2\right)}^{\left(0\right)}\times {I}_{3}}$ are core tensors. Note that the dimension of the third mode is retained without reduction. The reconstructed threeway array shown in Eq. 11 is the denoised grouping. Figure 4 illustrates the dimension reduction of GNDPCA with a threeway array.
After all the threeway arrays on the same level are denoised with GNDPCA, we place them back into their original positions. With the overlapping patches extracted, more than one patch estimation in the denoised threeway 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 ${\widehat{G}}_{0,0}$ of G _{0,0} is computed by an average of the patch estimates ${\widehat{\mathbf{Y}}}_{p}^{\left(0\right)}\in {R}^{{I}_{1}\times {I}_{2}},p=1,2,\dots ,P$ in the first stack. Here, ${\widehat{\mathbf{Y}}}_{p}^{\left(0\right)}\in {\widehat{\mathbf{Y}}}_{q}^{\left(0\right)},P\le Q$. We have
where E has the same size as ${\widehat{\mathbf{Y}}}_{p}^{\left(0\right)}$ and all elements of E are one, and x _{ R } is the coordinate iterated in patches. The denoised estimations for all stacks are denoted by ${\widehat{G}}_{a,b},0\le a\le A,0\le b\le 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 highfrequency 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 $\widehat{g}$ can be reconstructed from L _{ a } by,
and then,
In the reconstructed process, the lowfrequency information of G _{ a } is discarded and replaced by G _{ a + 1}. Finally, the despeckled image in logdomain is converted into spatial domain through exponential transformation. The complete MRGNDPCA method is shown as a diagrammatic sketch in Figure 2.
Experiments and results
The proposed MRGNDPCA 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 MRGNDPCA 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 lowpass filtering a complex Gaussian random field and taking the magnitude of the filtered output [27]. Lowpass filtering was performed by averaging the complex values of the size three sliding window. This shortterm 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 noisefree 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 MRGNDPCA method was compared with traditional despeckling methods (Median filter, Fourier ideal filtering, Fourier Butterworth filter, Wavelet filter) as well as stateoftheart algorithms (Bilateral filter [28], KSVD [12], Nonlocal 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 wellknown metrics as quantitative performance measures to assess quality image, including meansquare error (MSE), signaltonoise ratio (SNR), peak signaltonoise ratio (PSNR), and the structural similarity (SSIM). MSE evaluates the noise variance and is given by
where $\widehat{g}\left(x,y\right)$ 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 noisefree reference image as
PSNR is defined by
where N _{max} is the maximum fluctuation in the input image.
And SSIM is given by
where ${\mu}_{\widehat{g}\left(x,y\right)}$, μ _{ S(x,y)}, ${\sigma}_{\widehat{g}\left(x,y\right)}$, σ _{ S(x,y)}, and ${\sigma}_{\widehat{g}\left(x,y\right)S\left(x,y\right)}$ are the mean of $\widehat{g}\left(x,y\right)$, the mean of S(x, y), the variance of $\widehat{g}\left(x,y\right)$, the variance of S(x, y), and the covariance of $\widehat{g}\left(x,y\right)$ and S(x, y), respectively; C _{1} and C _{2} are small constants.
Multiple parameter selection
In MRGNDPCA, the parameters: number of modesubspace bases D _{ i(1)} × D _{ i(2)} , number of clusters (K), and the Gaussian scale (σ) are crucial for restoration results. In order to determine the best parameters of MRGNDPCA, we investigate the variation trends of MSE, SNR, PSNR and SSIM when D _{ i(1)}, D _{ i(2)} ∈ [1, 5], K = 50 ~ 90, and $\sigma =\sqrt{1/2},1,\sqrt{3/2}$. The comparison results are illustrated in Figures 6, 7, 8 and 9. It is clear in Figures 6, 7 and 8 that the lowest MSE, the highest SNR and PSNR can be obtained when D _{ i(1)} × D _{ i(2)} = 1 × 1, K = 80 and $\sigma =\sqrt{3/2}$, which are also considered as the best parameters in the following experiments, though SSIM is a little bit higher with σ = 1 than that with $\sigma =\sqrt{3/2}$ in Figure 9.
Comparison of MRPCA and MRGNDPCA
In order to illustrate the effectiveness of GNDPCA for image denoising, we substitute GNDPCA for PCA in the proposed MRGNDPCA framework, which is called MRPCA. The restoration results for low, middle, and highlevel noise are shown in Figure 10. All the parameters in MRPCA are set completely same as those in MRGNDPCA. The first row ((a), (c) and (e)) presents the restoration results obtained with MRGNDPCA, and the second row ((b), (d) and (f)) presents the restoration results obtained with MRPCA. The calculated MSE, SNR, PSNR and SSIM related to the original image are also given in the captions below the figure. We can see that MRGNDPCA gives better overall performance in terms of both noise suppression and detail preservation. However, MRPCA 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 and 3 clearly show that the proposed MRGNDPCA can obtain the lowest MSE, as well as the highest SNR and PSNR. Moreover, the SSIMs of the proposed method are nearest to one among all the denoising filters. Thus, the MRGNDPCA is an effective denoising method for reducing the noise as well as preserving the details. Comparisons of Figures 11, 12 and 13 show that the proposed MRGNDPCA can obtain more impressive denoising results for the image with high level noise. Although there are still some artifacts in the smoothed images with all the denoising methods, such as upper right object, the proposed MRGNDPCA 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 GNDPCA performance is an appropriate combination for ultrasound images with massive noise compare to the stateofart methods.
Phantom data validation
Phantom data called CIRS MODEL 057A are available at http://www.cirsinc.com/. The Bscan image of these phantom data was obtained from a commercial ultrasonic scanner with a 3.5 MHz phased array probe. The 470 × 620 regions of interest of the original and processed images are shown in Figures 14(a) and (b). The highlighted lines in Figures 14(a) and (b) indicate the 340th rows in the images. The pixel values in these rows from the original and denoised images are compared in Figures 14(c) and (d) respectively. Our proposed method suppresses speckle noise effectively. Most of the stray mottles in the original image are smoothed.
Clinical data verification
Clinical images of common carotid artery (CCA) are available for free on the Signal Processing Laboratory website. The Bscans were acquired with a Sonix OP ultrasound scanner with a different setup 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 MRGNDPCA. 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 MRGNDPCA 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 GNDPCA 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 GNDPCA 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 MRGNDPCA 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.
References
 1.
Tay PC, Garson CD, Acton ST, Hossack JA: Ultrasound despeckling for contrast enhancement. IEEE Trans Image Process 2010,19(7):1847–1860.
 2.
Adam D, BeilinNissan S, Friedman Z, Behar V: The combined effect of spatial compounding and nonlinear filtering on the speckle reduction in ultrasound images. Ultrasonics 2006,44(2):166–181. 10.1016/j.ultras.2005.10.003
 3.
Xavier FP, de Fontes GAB, Coupé P, Hellier P: Real time ultrasound image denoising. J RealTime Image Process 2011,6(1):15–22. 10.1007/s1155401001585
 4.
Kotropoulos C, Pitas I: Optimum nonlinear signal detection and estimation in the presence of ultrasonic speckle. Ultrason Imaging 1992,14(3):249–275. 10.1177/016173469201400303
 5.
Bernstein R: Adaptive nonlinear filters for simultaneous removal of different kind of noise in images. IEEE Trans Circuits Syst 1987,34(11):1275–1291. 10.1109/TCS.1987.1086066
 6.
Karaman M, Kutay MA, Bozdagi G: An adaptive speckle suppression filter for medical ultrasonic imaging. IEEE Trans Med Imaging 1995,14(2):283–292. 10.1109/42.387710
 7.
Yu Y, Acton ST: Speckle reducing anisotropic diffusion. IEEE Trans Med Imaging 2002,11(11):1260–1270.
 8.
Sanches JM, Marques JS: A MAP estimation algorithm with Gibbs prior using an IIR filter. energy minimization methods in computer vision and pattern recognition. Lect Notes Comput Sci 2003, 2683: 436–449. 10.1007/9783540450634_28
 9.
Portilla J, Strela V, Wainwright MJ, Simoncelli EP: Image denoising using scale mixtures of gaussians in the wavelet domain. IEEE Trans Image Process 2003,12(11):1338–1351. 10.1109/TIP.2003.818640
 10.
Buades A, Coll B, Morel JM: A review of image denoising algorithms, with a new one. Multiscale Model Simul 2005,4(2):490–530. 10.1137/040616024
 11.
Foi A, Boracchi G: Foveated selfsimilarity in nonlocal image filtering. In Proc. SPIE Electronic Imaging 2012. Burlingame (CA), USA: Human Vision and Electronic Imaging; 2012:8232–8291.
 12.
Aharon M, Elad M, Bruckstein AM: The KSVD: an algorithm for designing of overcomplete dictionaries for sparse representation. IEEE Trans Signal Process 2006,54(11):4311–4322.
 13.
Dabov K, Foi A, Katkovnik V, Egiazarian K: BM3D Image Denoising with Shape  Adaptive Principal Component Analysis. In Proc. Workshop on Signal Processing with Adaptive Sparse Structured Representations (SPARS'09). SaintMalo, France; 2009.
 14.
Zhang L, Dong W, Zhang D, Shi G: Twostage image denoising by principal component analysis with local pixel grouping. Pattern Recognit 2010, 43: 1531–1549. 10.1016/j.patcog.2009.09.023
 15.
Deledalle CA, Salmon J, Dalalyan A Proceedings of the British Machine Vision Conference. In Image denoising with patch based PCA: local versus global. Dundee, UK: BMVA Press; 2011:1–10.
 16.
Chen YW, Han X, Nozaki S: ICAdomain filtering for penumbral imaging. Rev Sci Instrum 2004, 75: 3977–3979. 10.1063/1.1787932
 17.
Han XH, Chen YW: A robust method based on ICA and mixture sparsity for edge detection in medical images for edge detection in medical images. Signal Image Video Process 2011,5(1):39–47. 10.1007/s1176000901405
 18.
Xu R, Chen YW: Generalized N dimensional principal component analysis (GNDPCA) and its application on construction of statistical appearance models for medical volumes with fewer samples. Neurocomputing 2009, 72: 10–12.
 19.
Jolliffe IT: Principal Component Analysis. SpringerVerlag; 1986:487. doi:10.1007/b98835
 20.
Kolda TG, Bader BW: Tensor decompositions and applications. SIAM Rev 2009,51(3):455–500. 10.1137/07070111X
 21.
Jain AK: Fundamental of Digital Image Processing. Englewood Cliffs, NJ: PrenticeHall; 1989.
 22.
Zong X1, Laine AF, Geiser EA: Speckle reduction and contrast enhancement of echocardiograms via multiscale nonlinear processing. IEEE Trans Med Imaging 1998,17(4):532–540. 10.1109/42.730398
 23.
Lowe DG: Distinctive image features from scaleinvariant keypoints. Int J Comput Vis 2004,60(2):91–110.
 24.
Kanungo T, Mount DM, Netanyahu NS, Piatko CD, Silverman R, Wu AY: An efficient kmeans clustering algorithm: analysis and implementation. IEEE Trans Pattern Anal Mach Intell 2002,24(7):881–892. 10.1109/TPAMI.2002.1017616
 25.
Jain AK: Data clustering: 50 years beyond Kmeans. Pattern Recognit Lett 2010,31(8):651–666. 10.1016/j.patrec.2009.09.011
 26.
Lu H, Plataniotis KN, Venetsanopoulos AN: MPCA: Multilinear Principal Component Analysis of tensor objects, IEEE Trans. Neural Netw 2008,19(1):18–39.
 27.
Piˇzurica A, Philips W, Lemahieu I, Acheroy M: A versatile wavelet domain noise filtration technique for medical imaging. IEEE Trans Med Imaging 2003,22(3):323–331. 10.1109/TMI.2003.809588
 28.
Tomasi C, Manduchi R: Bilateral filtering for gray and color images. Bombay: Sixth International Conference on Computer Vision; 1998:839–846.
 29.
Manj’on J, CarbonellCaballero J, Lull J, Garc’ıaMart’ı G, Mart’ıBonmat’ı L, Robles M: MRI denoising using nonlocal means. Med Image Anal 2008,12(4):514–523. 10.1016/j.media.2008.02.004
Acknowledgements
This work was supported by the National Basic Research Program of China (2010CB732505, 2013CB328806), the Key Projects in the National Science & Technology Pillar Program (2013BAI01B01), the National HiTech Research and Development Program (2013AA013703) and the National Science Foundation Program of China (61179020).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
ADN suggested the algorithm for images denoising, implemented it and drafted the manuscript. YJ conceived of the study, and participated in its design and coordination and helped to draft the manuscript. CY performed the acquisition of the ultrasound images and expressed opinions on the evaluation metric of the denoised results. CWJ and FJF participated in the design of the study, drew a part of figures and performed the image analysis. WYT helped to draft the manuscript and expressed opinions on the overall framework of the manuscript. All authors have read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Multiresolution
 Multilinear subspace learning
 Noise reduction