 Research
 Open Access
 Published:
Postreconstruction filtering of 3D PET images by using weighted higherorder singular value decomposition
BioMedical Engineering OnLinevolume 15, Article number: 102 (2016)
Abstract
Background
Positron emission tomography (PET) always suffers from high levels of noise due to the constraints of the injected dose and acquisition time, especially in the studies of dynamic PET imaging. To improve the quality of PET image, several approaches have been introduced to suppress noise. However, traditional filters often blur the image edges, or erase small detail, or rely on multiple parameters. In order to solve such problems, nonlocal denoising methods have been adapted to denoise PET images.
Methods
In this paper, we propose to use the weighted higherorder singular value decomposition for PET image denoising. We first modeled the noise in the PET image as Poisson distribution. Then, we transformed the noise to an additive Gaussian noise by use of the anscombe root transformation. Finally, we denoised the transformed image using the proposed higherorder singular value decomposition (HOSVD)based algorithms. The denoised results were compared with results from some general filters by performing physical phantom and mice studies.
Results
Compared to other commonly used filters, HOSVDbased denoising algorithms can preserve boundaries and quantitative accuracy better. The spatial resolution and the low activity features in PET image also can be preserved by use of HOSVDbased methods. Comparing with the standard HOSVDbased algorithm, the proposed weighted HOSVD algorithm can suppress the stairstep artifact, and the timeconsumption is about half of that needed by the Wieneraugmented HOSVD algorithm.
Conclusions
The proposed weighted HOSVD denoising algorithm can suppress noise while better preserving of boundary and quantity in PET images.
Background
Positron emission tomography (PET) is a noninvasive highspecificity imaging modality that is widely used for both clinical diagnosis and preclinical research because of its ability to reveal the metabolic evolution of molecular and biochemical processes in living organs. However, it suffers from relatively poor spatial resolution and high levels of noise due to the constraints of the injected dose and acquisition time, especially in the studies of dynamic PET imaging. High noise levels make the quantitative interpretation of PET images difficult; this is especially problematic when the image is to be used for early tumor detection when the tumor is small [1], myocardial perfusion [2], or hypometabolism of glucose in studies of early diagnosis of Alzheimer’s disease [3, 4]. In such applications, high noise level contained in PET images will cover small but important lesions and make the edges between organs blurred, further leading to mistakes in diagnosis and quantitation. Therefore, noise reduction is important for interpretation and computeraided analysis of PET images.
Several approaches have been introduced to suppress noise contained in PET images, during or after image reconstruction process. Statistical image reconstruction method such as maximum likelihood expectation maximization (MLEM) [5], its accelerated version ordered subset expectation maximization (OSEM) [6], or other modified algorithms [7–10] can model the noise in acquired data but suffer from increased noise as iterations increase. The regularization algorithm [11, 12], which reduces noise in the statistical reconstruction process by the way of adding a prior or penalty term to the reconstruction process, has promised a good denoising performance, but the determination of the appropriate regularization parameters is a challenging work, which is usually taskdependent. Alternatively, a smoothing filter [11, 13] after reconstruction, such as a Gaussian filter (GF), is widely used in current clinical applications. This kind of filter suppresses noise by averaging adjacent pixels by using a weight function of spatial distance between pixels [14], and has an inherent drawback of blurring edges, which may be clinically relevant. In [15], the authors used masked smoothing method for CT perfusion images denoising, but its performance for PET images has not been proved.
Numerous edgepreserving image filters, such as the bilateral filter (BF) [16, 17], anisotropic diffusion filters (ADFs) [18, 19], or priors [11, 20], which are able to detect and maintain edges while smoothing over noiseinduced fluctuations, have been introduced to denoise PET images. However, these approaches rely on multiple parameters specified by users according to the study to control the smoothness level and target intensity levels. Actually, it is impossible to preserve all of the interesting features simultaneously [21, 22]. Moreover, these filters usually erase small details, such as small tumors, which can lead to misdiagnosis. Algorithms that incorporate the anatomical information from CT or MRI into PET reconstruction [19, 23, 24] will improve the denoising performance, but the improvement of performance can only be achieved when the boundaries from anatomical image and PET image coincide with each other strongly; if not, these algorithms may introduce incorrect information leading to bias and artifacts [25]. Transformbased filters, such as waveletbased filters [26–28] and singular value decomposition [29], have also been proposed to denoise PET images by exploiting the uncorrelated properties of noise in the transform domain.
Recently, nonlocal denoising methods, such as the nonlocal mean (NLM) filters [30] and blockmatching 3D (BM3D) filter [31], have been adapted to denoise medical images [32–35], which exploit the similarity and sparsity among small patches to preserve details while removing noise, and achieve outstanding performance. In [21], the authors extended the application of the NLM filter to the 3D PET image and modified it by incorporating anatomical information from the CT image to constrain the similarity measure within a subset of voxels. This algorithm can produce improved lesion contrast and signaltonoise ratio (SNR) when compared with other traditional PET image denoising methods.
The higherorder singular value decomposition (HOSVD) [36] is another kind of nonlocal denoising method, which is a development on highorder matrix or tensor of the singular value decomposition (SVD) of the twodimension matrix. It offers a simple method for handling sparsity among similar patches by grouping them into a highorder matrix, then denoising the image in a transformthresholdinverse transform fashion. The transform bases are learned from the original image and are more adaptable to the image content and may achieve a much sparser representation than fixed bases used in the BM3D method.
The HOSVDbased method has been extended into the applications of MR image denoising, and demonstrated stateoftheart performance [37]. In this study, we aim to investigate and extend the application of HOSVD to denoise 3D PET images. We examined its ability to preserve highintensity features in denoised images, and the ability to detect potentially low remnant activity by the percentage recovered signal (PRS) in lesion volume of interest (VOI) and the contrasttonoise ratio (CNR) in the denoised images. Results were compared to those obtained after GF, BF, and ADF filtering of the same data.
Methods
To extend HOSVD to PET image denoising and propose our HOSVDbased 3D PET images denoising method, we will first give a brief introduction about HOSVD in image denoising. Then, we will introduce the Anscombe variancestabilizing transformation and its exact unbiased inverse. Finally, we will describe our method and its implementation.
Nonlocal HOSVD denoising for 2D images
Similar to most denoising algorithms, the HOSVD denoising method is originally designed for the additive zero mean Gaussian noise. The main assumption in the denoising process is that the observed noisy image Y from a noise free image X with additive Gaussian noise N can be modeled as \(Y = X+N\), where N has a known variance \({\sigma }^2\) and zero mean.
To find a good estimate \(\tilde{X}\) of X from Y, the HOSVD denoising method first groups the similar patches for each reference image patch into a 3D stack. Then, the HOSVD transformation of the stack is performed to obtain its representation by the HOSVD bases and coefficients. Based on the uncorrelated properties of noise in the HOSVD transform domain, the coefficients are truncated by a hard thresholding operator, then one denoised estimate of this stack is obtained by the inverse HOSVD transform. Repeating this operation for each reference patch results in multiple estimates for each pixel. The final denoised image is obtained by averaging these estimates. The details of the standard HOSVD denoising process are described in the following.
Given a \(p\,{\times }\,p\) reference patch \(P_{ref}\) in the original noisy 2D image, K similar patches are found including \(P_{ref}\) and stacked into a 3D matrix \(Z\,{\in}\,R^{p\,{\times }\,p{\times }\,K}\). Then, the HOSVD of Z is performed and can be formulated as [38]
where \(U^{(1)},U^{(2)}\,{\in}\,R^{p\,{\times }\,p}\), and \(U^{(3)}\,{\in}\,R^{K\,{\times }\,K}\) are orthonormal unitary matrices. They are computed from the SVD of the unfolding matrices of Z along three directions respectively [38]. \({\times }_n\) stands for the nmode tensor product. S is a 3D coefficient matrix of size \(p\,{\times}\,p\,{\times}\,K\). Eq. (1) indicates that the HOSVD can exploit signal sparsity across each dimension.
The coefficient matrix S can be computed as:
where T represents the transpose of a matrix. Then, stack Z can be filtered by a hardthreshold operation for S, and the threshold \(\tau\) can be determined as follows [39]:
Thus, the truncated coefficient matrix \(\tilde{S}\) with entries \(\tilde{s}_i\) can be determined as follows:
where the threshold \(\tau\) is optimal from a statistical risk viewpoint for any orthonormal basis and for \(N{\sim }\mathcal {N}(0,\sigma ^2)\). This truncation of the HOSVD coefficients is based on the assumption that the underlying noisefree image can be sparsely represented by the HOSVD bases. Using (1), we can obtain an estimate \(\tilde{Z}\) of Z by following:
Repeating the above process for each reference patch in sliding window fashion, multiple estimates at each pixel will be obtained. The final filtered image is produced by averaging these estimates at each pixel.
To improve the denoising performance even further, an additional Wiener filter step after the standard HOSVD denoising process can be performed. Let \(\tilde{Z}\) be a stack of similar patches from the standard HOSVD filtered image, which corresponds to Z from the original noisy image, then the coefficients of \(\tilde{Z}\) and Z in the HOSVD bases of \(\tilde{Z}\) are denoted as \(\tilde{s}\) and s respectively. Then, the filtered coefficients of Z, denoted as \(\bar{s}\), are computed as follows:
In this paper, the standard HOSVD and its Wiener filteraugmented version are referred to as HOSVDS and HOSVDW respectively.
HOSVDbased denoising for 2D PET images
The abovementioned HOSVD denoising for 2D images cannot be straightly extended to denoise PET images, as (1) for nonlinear reconstruction algorithms, usually, the exact noise model in reconstructed PET images is unknown, even though raw PET data follow the Poisson distribution, and (2) the intensity value at each voxel in PET images is the representation of the activity concentration (Bq/ml) in the corresponding region of the space, which is low, so the assumption of an imageindependent noise model does not hold.
In fact, for images with low photon counts, such as a PET image or SPECT image, the pixel or voxel intensity can be modeled approximatively with the Poisson distribution, which follows the observations that the noise variation in PET images are multiplicative, asymmetric and varying in the image. In this paper, we approximate noise in reconstructed PET images with the Poisson model. Then, with the Anscombe root transformation [40] and its inverse, the HOSVDbased denoising algorithms can be used to denoise PET images. Here, we employ the closedform approximation of the exact unbiased inverse of the Anscombe variancestabilizing transformation [41] to alleviate the biased estimates caused by the nonlinearity of the Anscombe transformation. The Anscombe transformation and its inverse are denoted as \(AVST\) and \(IAVST\), respectively, which are formulated as follows:
and
where I represents the noisy PET image, and Z is the denoised data in the HOSVD transform domain.
The HOSVDbased denoising algorithms can then be extended to the denoising application of PET images by using the following formula:
where \({Den}_{\scriptscriptstyle {HOSVD}}(\cdot )\) represents the HOSVD denoising manipulation, I and \(\tilde{I}\) are the noisy and denoised PET images, respectively.
Weighted HOSVD denoising scheme for 3D PET images
The hardthreshold operation used in HOSVDS may result quite different between reference and selected similar patches in structure; this is especially true when the exact structure is unavailable in the underlying clean images, as with PET image denoising. In such applications, straight use of HOSVDS may lead to stairstep artifact in the denoised image, as shown in Fig. 1. Although the HOSVDW can modify this disadvantage, its use also increases time consumption dramatically to about twice the time required for HOSVDS.
To improve the denoising performance of the HOSVD while keeping the time consumption from dramatically increasing, we assign different weights to different selected similar patches, depending on their level of similarity with the reference patch, in structure and statistics. This is different to the HOSVDS method, which assign identical weight to each selected similar patch. The similarity level in statistics is measured by the pvalue output of the KolmogorovSmirnov (KS) test between reference patch and the selected similar patch.
For the calculation of weights, we first consider the patch similarity measure used in the HOSVDS. The distance between the threedimension (3D) reference patch \(P_{ref}\) and its ith similar patch \(P_i\) can be formulated as:
Note that the reference patch \(P_{ref}\) and a patch P within the search window is similar only if the Euclidean distance \(d(P_{ref},P)\) between them meets the following rule [36]:
Based on this, we calculate the similarity weight between \(P_i\) and \(P_{ref}\) as follows:
which means that patches with high similarity to the reference patch will be assigned a larger weight.
To alleviate the disadvantage caused by the original similarity measure, we used a hypothesis test (the onesided KolmogorovSmirnov test) to check how well the values in \(P_{ref}P_i\) conform to \(\mathcal {N}(0,2\sigma ^2)\), and the pvalue output by the KS tests are used to modify the weighting factor by a multiplication factor. We denote the p value output by the KS test as \(p_{KS}(P_{ref},P_i)\), then the weight of patches that are similar to the reference patch \(P_{ref}\) can be formulated as follows:
Thus, the ith 3D patch that is similar to \(P_{ref}\) can be formulated as follows:
Then, the clustering 4D stack can be formulated as:
Note that \(C(P_{ref})\) is a matrix of order four. The denoised 4D stack is formulated as:
The estimation \(\breve{P}_i\) of \(P_i\) can then be obtained as follows:
Repeating the above process for each reference patch and averaging the estimates at each voxel, we can obtain the denoised image in the Anscombe transform domain. The final filtered PET image can then be produced by performing the \(IAVST (\cdot )\) operation. For clarity, the proposed weighted HOSVD denoising algorithm is referred to as HOSVDKSW, and the flow diagram is shown in Fig. 2.
HOSVDKSW not only preserves all the advantages of the HOSVDS method, but suppresses the stairstep artifact caused by the HOSVDS with similar time consumption.
Practical implementation of the HOSVDbased denoising algorithm
In this work, three HOSVDbased methods are applied to denoise 3D PET images. The most timeconsuming parts of HOSVDbased algorithms are 3D patch grouping and HOSVD transform. In the HOSVDS and HOSVDW implementation, the Knearestneighbors of the reference patch are found based on (11) and are then grouped into a stack. For HOSVDKSW, the patches that correspond to the first Kmaximumweight are selected as the Knearestneighbors of the reference patch (including the reference patch itself) and similar patches obtained by using (11).
To improve the computation efficiency while preserving the denoising efficacy of algorithms, we restricted the size of the searching space of similar patches and the sliding step of the reference patch. On one hand, bigger search windows do not always lead to better results, and a further increase of the search window will increase time consumption, especially in 3D PET image denoising, where similar structures of a special reference patch are usually concentrated in a restricted area close to it. On the other hand, the algorithm execution time is linearly corrected with the number of reference patches selected in the implementation.
In our implementation of the HOSVDbased algorithms, the size of the search window is set to 20 voxels and the sliding step is set to \(p1\) at each image dimension, where p is the size of the reference patch in each dimension.
Experiments and results
Phantom study
To assess various effects in image quality resulting from different denoising methods, two physical phantom studies were performed, one for the miniDerenzo phantom and another for the small animal image quality phantom which is described by NEMA NU 42008 [42].
The NEMA small animal image quality phantom (homemade), created from polymethylmethacrylate, has three main parts, as shown in Fig. 3a. The first is a homogeneous cavity with a 30 mm diameter and 30 mm length, filled with radioactivity (hot region) to measure the uniformity. The second part consists of two chambers housed in the cavity and filled with water and air (cold regions, 15 mm in length, wall thickness and outer diameter are 1 and 10 mm respectively) to estimate the spillover ratio. The third part of the phantom contains five rods which are drilled through a plastic region and have diameters of 1–5 mm respectively, and are used to measure the noise and recovery coefficients.
The miniDerenzo phantom with rod diameters ranging from 0.7 to 2.4 mm was used to compare the spatial resolution characteristics of the different methods. Fig. 3b is the schematic diagram of the miniDerenzo phantom.
The small animal image quality phantom was filled with 2.74 MBq of 18FFDG, and the miniDerenzo phantom had 1.11 MBq. Data were acquired for 15 min for both phantoms using a GENISYS\(^4\) PET scanner (Sofie Biosciences Inc., USA). The images were reconstructed on the GENISYS\(^4\) Acquisition Engine with 3D OSEM algorithm (without PSF) into 96 × 96 × 208 voxels, where the voxel dimensions were 0.457 × 0.457 × 0.457 mm.
Nude mice 18FFDG study
To compare the different denoising approaches in terms of their effect on preclinical or clinical quantitative studies, we performed an 18FFDG PET study of nude mice with 4T1 tumors. A 2.24 MBq injection of 18FFDG mixed with isotonic saline to a total volume of 150 \(\upmu\)L was administered into nude mice through a tail vein with the mouse immobilized in a restrainer. Image acquisition started 45 min after injection and was performed on the GENISYS\(^4\) PET scanner without attenuation correction.
Quantitative evaluations
Uniformity
To quantitatively evaluate and compare the effect of different denoising algorithms on the uniformity of the PET image, a 22.5 mm diameter (75 % of the active diameter) by 10 mm long cylindrical volume of interest (VOI) was drawn over the center of the uniform region of the NEMA small animal image quality phantom. The average activity concentration in this VOI, and its percentage standard deviation (%STD) were calculated. The percentage standard deviation was calculated as follows:
where
is the standard derivation of activity concentration in a uniform region VOI, and \({Act}_i\) is the activity concentration in the ith voxel of the VOI, and n represents the number of voxels that are contained in uniform region VOI. \(STD_{uniform}\) and \(Mean_{uniform}\) are the standard deviation (STD) and mean in the uniform region VOI, respectively.
Recovery coefficients
A single image slice with lower noise characteristics was obtained by averaging image slices that cover the central 10 mm length of the rods. Circular region of interests (ROIs) were taken in this image, around each rod with diameters twice the physical diameter of the rods [42]. Transverse image pixel coordinates of the locations with the maximum voxel value within each of these ROIs were used to create line profiles along the rod in the axial direction. Then, the mean and standard deviation of the recovery coefficient (RC) for the rods were calculated as follows:
where \(Mean_{line}\) represents the mean of the activity concentration for each line profile. The percentage standard deviation of the recovery coefficient for each rod was calculated as follows:
Spillover ratio
For the calculation of the spillover ratio of water and airfilled regions, two cylindrical VOIs with a diameter of 4 mm (half the physical diameter of the cylinders) and encompassing the central 7.5 mm in length were drawn in the water and airfilled cylindrical inserts respectively. Then, the spillover ratio (SOR) was calculated as follows:
where \(Mean_{VOI}\) was the mean of the activity concentration in the water or airfilled region VOI. The percentage standard deviation of SOR was calculated as follows:
Spatial resolution
The influence of different algorithms to the spatial resolution of the PET image was evaluated visually by the miniDerenzo phantom study.
Percentage recovery signal and contrasttonoise ratio
In the nude mice study, we compared the ability to preserve highintensity features of different denoising approaches by two metrics. The first metric is the PRS in a VOI in the denoised image relative to noisy image, and the second is the CNR. These two metrics are important in clinical diagnosis and treatment monitoring. In comparison to PRS, which fails to capture the noise characteristics of the noisy images, CNR is more holistic, especially in studies where the goal is to detect potentially low remnant activity, the CNR can be quite critical.
Two VOIs were drawn for the calculation of PSR and CNR, a hot VOI in the lesion region and a lowintensity VOI that has a relatively uniform uptake in the background region (as shown in Fig. 12a). The PRS and CNR were calculated as follows:
and
where \(M_{les}^{de}\) and \(M_{les}^{un}\) denote the mean intensity in the lesion VOI in the denoised images and noisy images respectively. \(M_{les}\) denotes the mean intensity in the lesion VOI in denoised or noisy images, and \(M_{back}\) and \(STD_{back}\) represent the mean and standard deviation of intensity in the low intensity (background) VOI in the same images.
Filtering parameters
In this work, 3D GF, BF, and ADF were applied. The only parameter of GF is the standard deviation of the Gaussian kernel, which is denoted as \(\sigma _d\) and used to control the diffusion rate of the Gaussian kernel. For BF, there are two parameters to be determined to control the smoothing degree of the images, namely the geometric spread \(\sigma _d\) and photometric spread \(\sigma _r\), which determine the geometric closeness and photometric similarity of the neighboring voxels to be used in the smoothing process, in domain and range respectively. There are many flow functions which can be used in ADF, however, considering that one purpose of PET image denoising is to preserve a prominent signal in a homogeneous area, such as a focal lesion within a single anatomical region, we chose the exponential flow function in this work. The smoothing performance of ADF was controlled by two parameters, the diffusion rate \(\lambda\), and the edge preserving parameter \(\kappa\). The performance of the HOSVDbased algorithms depends on the setting of two patchrelated parameters, patch size p and the number K of similar patches in a group; however, in a rational range, K would not significantly affect the denoising performance. Therefore p is the only free parameter which influences the accuracy in HOSVDbased algorithms.
In order to compare the performance of HOSVDKSW with GF, BF, and ADF for PET image denoising, filtering parameters that make all these algorithms balanced should be determined; these selected parameters result similar standard deviation in the homogeneous regions in denoised images. In the phantom study, we chose this region as the uniform VOI, and in the mice study, the background VOI was used. The comparison between HOSVDbased algorithms was performed by using the same patch size p. For GF, the kernel size is fixed to five voxels. We fixed \(\sigma _d\) at 4.5 in BF. For ADF, keeping \(\lambda = 0.15\) (\(\lambda \in [0,0.25]\) should be sufficiently small to make the diffusion process stable [18]) unchanged. The residual parameters were swept over the range listed in Table 1. Note that the parameter \(\kappa\) in ADF was represented as the percentage of standard deviation inside the selected homogeneous region in this work.
Figure 4 shows the variation of the standard deviation in uniformity region VOI when different filtering parameters setting were used in the phantom study. The parameters corresponding to the standard deviation value circled in red dashed circles were chosen in the phantom study. For GF, \(\sigma _d\,=\,1\) (1.07 mm FWHM) was used. For BF, \(\sigma _d\,=\,4.5\) and \(\sigma _r\,=\,25.5\). \(\lambda\,=\,0.15\) and \(\kappa\,=\,75\,\%\) for ADF. For HOSVDS, HOSVDW, and HOSVDKSW, patch size was set to \(p\,=\,7\).
For the nude mice study, the filtering parameters were \(\sigma _d\,=\,1\) for GF, \(\sigma _d\,=\,3.0\) and \(\sigma _r\,=\,10\) for BF, \(\lambda\,=\,0.15\) and \(\kappa\,=\,80\,\%\) for ADF, and \(p\,=\,9\) for HOSVDS, HOSVDW, and HOSVDKSW.
Results
Image quality
Figure 5 shows the mean (Fig. 5a) and percentage standard deviation (Fig. 5b) of the noisy and denoised image in the uniform VOI. From this result, all of these filters can preserve the mean and improve the smoothness in the uniform region. However, HOSVDKSW yielded a similar smoothing result with GF, BF, and ADF, which were smoother than HOSVDS and HOSVDW. As was expected, GF yielded an oversmoothed boundary in the uniform region while BF and ADF preserved boundary. A similar edgepreserving result also could be observed in HOSVDbased methods. This can be seen in Fig. 6.
Figure 7 shows how the recovery coefficients (RC) of rods varies with diameters in the noisy and denoised image under different algorithms. From Fig. 7a, the GF, BF, and ADF deteriorated the RC over all of the rods, while the HOSVDbased algorithms could preserve the RC, and the HOSVDKSW improved RC slightly even with a smaller rod diameter. Figure 7b shows that HOSVDS is less stable when compared with HOSVDKSW, even though the Wiener filter step in HOSVDW would improve this shortage in HOSVDS.
Figure 8 shows the results of the spillover ratio and the corresponding percentage standard deviation for the cavities filled with water and air. As shown in Fig. 8a, the HOSVDKSW decreased SORs both in water and airfilled regions while other compared methods enlarged it besides the GF. However, compared with HOSVDKSW, the improvement of GF to SOR in water and airfilled regions was negligible. Again, the percentage standard deviation (%STD) of SOR in water and airfilled regions resulting from the HOSVDKSW are comparable to those from GF and BF, and the HOSVDS and HOSVDW resulted in a higher %STD than other algorithms although both of them reduced it when compare to the noisy image, this was shown in Fig. 8b. Figure 9 shows the profiles of the dotted lines drawn in the center of the water/air chambers region, a reduction of the normalized counts in cold chambers could be observed in HOSVDKSW.
MiniDerenzo phantom
The results of the miniDerenzo phantom study are shown in Fig. 10. The noisy image is shown in Fig. 10a. The denoised image using GF with \(\sigma _d=1\) is shown in Fig. 10b. As expected, the GF smoothed out signal and noise simultaneously. The BF denoised image demonstrated better preserved hot boundaries than the GF, but the cold boundaries were smoothed out, as denoted by the red arrow in Fig. 10c. ADF yielded a result between GF and BF as shown in Fig. 10d. HOSVDbased algorithms demonstrated betterpreserved results of cold boundaries and hot regions. There is no evident difference in the denoised images resulting from HOSVDS, HOSVDW, and HOSVDKSW. This may have benefit from the accurate structural information existing in the miniDerenzo phantom, as the simple similarity measure used in HOSVDS is just enough for similar patch selection for a certain reference patch.
As shown in Fig. 10, the activityfilled area with the diameter of 1.6 mm can be distinguished clearly in all denoised images. However, the area with a diameter of 1.2 mm was smoothed in GF, BF, and ADF, while the HOSVDbased algorithms demonstrated betterpreserved details in this area. It is clear that the HOSVDbased algorithms can preserve the spatial resolution of the PET system well.
Nude mice study result
Figure 11 presents coronal slices through lesions in the nude mice 18FFDG PET study. PET image shown in (a) is a noisy reconstructed image, (b), (c), (d), (e), (f), and (g) are the images that are denoised by GF, BF, ADF, HOSVDS, HOSVDW, and HOSVDKSW, respectively. Lesions are denoted by yellow arrows in these images. The lesion is ellipsoidal with its length in three axis are 1.1, 0.9 and 0.6 cm respectively, and the volume is 2.5 mL. It can be observed that GF suppressed noise as well as the lesion and organ boundaries. BF suppressed noise to a similar extent as GF in the background region, and yielded sharper boundaries of the lesion and organs, as shown in Fig. 11c. ADF preserved the high contrast objects, but low contrast details such as the lesion were smoothed, and there were stairstep artifacts in the boundaries of high contrast objects as denoted by red arrows shown in Fig. 11d. All of the HOSVDbased algorithms suppressed noise effectively, however, HOSVDS yielded stairstep artifacts, as shown in Fig. 11e. Both HOSVDW and HOSVDKSW preserved the boundaries of the lesion and other high contrast objects, as can be seen in Fig. 11f, g.
Quantitative results from the nude mice 18FFDG PET study are shown in Fig. 12. VOIs in the lesion and background regions used for computing the CNR and PSR are delineated in Fig. 12a. It can be seen that all compared algorithms improved the CNR (Fig. 12b) in the lesion region in comparison to the noisy image. However, GF, BF, and ADF yielded worse results than HOSVDbased algorithms. GF improved CNR by 249 % (CNR = 18) compared to the noisy image (CNR = 5), and introduced 9 % (PSR = 91 %) reduction in average normalized counts in the lesion. BF and ADF improved CNR by 208 and 200 %, and yielded 8 % (PRS = 92 %) and 9 % (PNS = 91 %) reduction in the average normalized counts in the same region, respectively. These demonstrated that traditional edge preserving filters could not achieve better quantitative performance than GF when their smoothness of the background region was similar. However, HOSVDbased algorithms introduced a better CNR and yielded better PRS than other compared algorithms, among which HOSVDW yielded the best performance, with a 387 % increase in CNR, while maintaining similar average normalized counts (PSR = 96 %) in the lesion when compare to the other two HOSVDbased methods. Compared to HOSVDW, HOSVDKSW yielded the second best CNR and PRS with CNR = 24 and PSR = 96 %, but its time consumption was dramatically lower than for the HOSVDW.
Discussion
In this study, we investigated the use of a novel denoising algorithm based on highorder singular value decomposition in PET image denoising. This algorithm performs image denoising by exploiting the sparsity between similar patches by using the HOSVD transform, and offers a simple and elegant method for handling sparsity among similar patches. The HOSVD bases are learned from the image and thus are more adaptable to the image content and may achieve a much sparser representation than fixed bases such as wavelet and discrete cosine bases.
To evaluate the performance of HOSVDbased denoising algorithms in PET image filtering, we performed three 18FFDG PET studies with: a NEMA small animal image quality phantom, a miniDerenzo phantom, and a nude mice with 4T1 tumors. The noisy reconstructed PET images were filtered by GF, BF, ADF, and HOSVDbased denoising methods. Our results demonstrated that HOSVDbased methods yielded consistently superior performance compared with GF, BF, and ADF, in terms of visual quality and quantitative metrics. HOSVDW and HOSVDKSW recovered the small hot region better, while preserving the spatial resolution and reducing the SOR of the background regions. On condition that the smoothness of the background region is similar with each other, both methods suppressed the noise in soft tissues effectively and preserved the organ boundaries.
The straight usage of HOSVDS might lead to stairstep artifacts in denoised PET images. This is because the patch similarity measure may collect patch pairs with quite different structures in the Anscombe transformed image with identical weight. If these patches are grouped into a 4D stack then the HOSVD bases that learned from them are not accurate, and lead to inappropriate information spreading into the final denoised PET image. While the Wiener filter step used in HOSVDW can modify this shortcoming, the time consumption will dramatically increase to about twice as long as the HOSVDS. Our results demonstrated that the proposed method HOSVDKSW yielded similar performance compared with HOSVDW, and the timeconsumption was similar with HOSVDS.
The good denoising performance of the HOSVDKSW is attributed to two main features. The first feature is that HOSVDKSW groups the similar patches into a stack with a weighted manner, where the weight is calculated according to the similarity between reference patch and special patch in domain and statistics. This will make the HOSVD bases more accurate than those in HOSVDS, and thus remove the stairstep artifacts existing in HOSVDS denoised images. The second feature is that the similarity weight of similar patches is calculated by considering two factors. The first factor is the Euclidean distance between a similar patch and reference patch from the similarity measure process, and the second factor is the p value output of the KS test of the difference patch between the reference patch and the similar patch. This twopart criterion enables the high accuracy of the HOSVD bases.
The main shortcoming of the proposed HOSVDKSW method for 3D PET image denoising is that we considered the noise in the reconstructed PET images as the Poisson distribution. It was not accurate enough, as well known, and a more accurate noise model and corresponding variancestabilizing transformation should be used to improve the reliability of the results.
Conclusions
In summary, our study demonstrated that the HOSVDW and the proposed HOSVDKSW methods yield improved image quality while preserving the accuracy of lesion quantification. Considering that the timeconsumption of the HOSVDKSW is about half of the HOSVDW, we believe that the HOSVDKSW is more practical than the HOSVDW at a very low performance compromise.
Abbreviations
 PET:

positron emission tomography
 HOSVD:

higherorder singular value decomposition
 AVST:

anscombe variancestabilizing transformation
 MLEM:

maximum likelihood expectation maximization
 OSEM:

ordered subset expectation maximization
 GF:

Gaussian filter
 BF:

bilateral filter
 ADF:

anisotropic diffusion filters
 NLM:

nonlocal mean
 BM3D:

blockmatching three dimension filter
 SNR:

signaltonoise ratio
 SVD:

singular value decomposition
 PRS:

percentage recovered signal
 VOI:

volume of interest
 CNR:

contrasttonoise ratio
 HOSVDS:

standard HOSVD
 HOSVDW:

wiener filteraugmented HOSVD
 KS test:

KolmogorovSmirnov test
 HOSVDKSW:

weighted HOSVD
 PSF:

point spread function
 STD:

standard deviation
 ROI:

region of interest
 RC:

recovery coefficient
 SOR:

spillover ratio
References
 1.
Zhao C, Chen Z, Ye X, Zhang Y, Zhan H, Aburano T, Tian M, Zhang H. Imaging a pancreatic carcinoma xenograft model with 11Cacetate: a comparison study with 18FFDG. Nucl Med Commun. 2009;30(12):971–7.
 2.
Chen K, Lawson M, Reiman E, Cooper A, Feng D, Huang SC, Bandy D, Ho D, Yun LS, Palant A. Generalized linear least squares method for fast generation of myocardial blood flow parametric images with N13 ammonia PET. IEEE Trans Med Imaging. 1998;17(2):236–43.
 3.
Kwan RKS, Evans AC, Pike GB. MRI simulationbased evaluation of imageprocessing and classification methods. IEEE Trans Med Imaging. 1999;18(11):1085–97.
 4.
Silverman DH, Small GW, Chang CY, Lu CS, Kung De Aburto MA, Chen W, Czernin J, Rapoport SI, Pietrini P, Alexander GE, Schapiro MB, Jagust WJ, Hoffman JM, WelshBohmer KA, Alavi A, Clark CM, Salmon E, de Leon MJ, Mielke R, Cummings JL, Kowell AP, Gambhir SS, Hoh CK, Phelps ME. Positron emission tomography in evaluation of dementia: regional brain metabolism and longterm outcome. JAMA. 2001;17(286):2120–7.
 5.
Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging. 1982;1(2):113–22.
 6.
Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging. 1994;13(4):601–9.
 7.
Ma J, Hudson M. Blockiterative fisher scoring algorithms for maximum penalized likelihood image reconstruction in emission tomography. IEEE Trans Med Imaging. 2008;27(8):1130–42.
 8.
Sitek A. Reconstruction of emission tomography data using origin ensembles. IEEE Trans Med Imaging. 2011;30(4):946–56.
 9.
SzirmayKalos L, Magdics M, Toth B, Bukki T. Averaging and metropolis iterations for positron emission tomography. IEEE Trans Med Imaging. 2013;32(3):589–600.
 10.
Vaissier PEB, Goorden MC, Taylor AB, Beekman FJ. Fast countregulated OSEM reconstruction with adaptive resolution recovery. IEEE Trans Med Imaging. 2013;32(12):2250–61.
 11.
Nuyts J, Fessler JA. A penalizedlikelihood image reconstruction method for emission tomography, compared to postsmoothed maximumlikelihood with matched spatial resolution. IEEE Trans Med Imaging. 2003;22(9):1042–52.
 12.
Chen J, Tu K, Chen T, Lu H, Liu R, Chou K, Chen C. Crossreference maximum likelihood reconstruction for positron emission tomography with empiric studies. Eur J Nucl Med. 2001;28(8S):1073.
 13.
Nuyts J, Baete K, Beque D, Dupont P. Comparison between MAP and postprocessed ML for image reconstruction in emission tomography when anatomical knowledge is available. IEEE Trans Med Imaging. 2005;24(5):667–75.
 14.
Qi J, Leahy RM. Iterative reconstruction techniques in emission computed tomography. Phys Med Biol. 2006;51(15):541–78.
 15.
Wack DS, Snyder KV, Seals KF, Siddiqui AH. Masked smoothing using separable kernels for ct perfusion images. BMC Med Imaging. 2014;14:25.
 16.
Tomasi C, Manduchi R. Bilateral filtering for gray and color images. In: Alley G, editors. Sixth international conference on computer vision, 1998. Bombay: IEEE; 1998. p. 839–46.
 17.
Hofheinz F, Langner J, BeuthienBaumann B, Oehme L, Steinbach J, Kotzerke J, van den Hoff J. Suitability of bilateral filtering for edgepreserving noise reduction in PET. EJNMMI Res. 2011;1(1):23.
 18.
Perona P, Malik J. Scalespace and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell. 1990;12(7):629–39.
 19.
Chan C, Fulton R, Feng DD, Meikle S. Regularized image reconstruction with an anatomically adaptive prior for positron emission tomography. Phys Med Biol. 2009;54(24):7379–400.
 20.
Mumcuoğlu EU, Leahy RM, Cherry SR. Bayesian reconstruction of PET images: methodology and performance analysis. Phys Med Biol. 1996;41(9):1777–807.
 21.
Chan C, Fulton R, Barnett R, Feng DD, Meikle S. Postreconstruction nonlocal means filtering of wholebody PET with an anatomical prior. IEEE Trans Med Imaging. 2014;33(3):636–50.
 22.
Nuyts J, Michel C, Brepoels L, De Ceuninck L, Deroose C, Goffin K, Mottaghy FM, Stroobants S, Van Riet J, Verscuren R. Performance of MAP reconstruction for hot lesion detection in wholebody PET/CT: an evaluation with human and numerical observers. IEEE Trans Med Imaging. 2009;28(1):67–73.
 23.
Vunckx K, Atre A, Baete K, Reilhac A, Deroose CM, Van Laere K, Nuyts J. Evaluation of three MRIbased anatomical priors for quantitative PET brain imaging. IEEE Trans Med Imaging. 2012;31(3):599–612.
 24.
Somayajula S, Panagiotou C, Rangarajan A, Li Q, Arridge SR, Leahy RM. PET image reconstruction using information theoretic anatomical priors. IEEE Trans Med Imaging. 2011;30(3):537–49.
 25.
Comtat C, Kinahan PE, Fessler JA, Beyer T, Townsend DW, Defrise M, Michel C. Clinically feasible reconstruction of 3D wholebody PET/CT data using blurred anatomical labels. Phys Med Biol. 2002;47(1):1–20.
 26.
Stefan W, Chen K, Guo H, Renaut RA, Roudenko S. Waveletbased denoising of positron emission tomography scans. J Sci Comput. 2012;50(3):665–77.
 27.
Shidahara M, Ikoma Y, Kershaw J, Kimura Y, Naganawa M, Watabe H. PET kinetic analysis: wavelet denoising of dynamic PET data with application to parametric imaging. Ann Nucl Med. 2007;21(7):379–86.
 28.
Le Pogam A, Hanzouli H, Hatt M, Cheze Le Rest C, Visvikis D. Denoising of PET images by combining wavelets and curvelets for improved preservation of resolution and quantitation. Med Image Anal. 2013;17(8):877–91.
 29.
Bagci U, Mollura DJ. Denoising PET images using singular value thresholding and stein’s unbiased risk estimate. Med Image Comput Comput Assist Interv. 2013;16(Pt 3):115–22.
 30.
Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. Multiscale Model Simul. 2005;4(2):490–530.
 31.
Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3D transformdomain collaborative filtering. IEEE Trans Image Process. 2007;16(8):2080–95.
 32.
Manjón JV, Coupé P, Buades A, Louis Collins D, Robles M. New methods for MRI denoising based on sparseness and selfsimilarity. Med Image Anal. 2012;16(1):18–27.
 33.
Manjón JV, Coupé P, MartíBonmatí L, Collins DL, Robles M. Adaptive nonlocal means denoising of MR images with spatially varying noise levels. J Magn Reson Imaging. 2010;31(1):192–203.
 34.
Aksam Iftikhar M, Jalil A, Rathore S, Hussain M. Robust brain MRI denoising and segmentation using enhanced nonlocal means algorithm. Int J Imaging Syst Technol. 2014;24(1):52–66.
 35.
Maggioni M, Katkovnik V, Egiazarian K, Foi A. Nonlocal transformdomain filter for volumetric data denoising and reconstruction. IEEE Trans Image Process. 2013;22(1):119–33.
 36.
Rajwade A, Rangarajan A, Banerjee A. Image denoising using the higher order singular value decomposition. IEEE Trans Pattern Anal Mach Intell. 2013;35(4):849–62.
 37.
Zhang X, Xu Z, Jia N, Yang W, Feng Q, Chen W, Feng Y. Denoising of 3D magnetic resonance images by using higherorder singular value decomposition. Med Image Anal. 2015;19(1):75–86.
 38.
Lathauwer LD. Signal processing based on multilinear algebra. Leuven: Katholike Universiteit Leuven; 1997 (PhD thesis).
 39.
Donoho DL, Johnstone IM. Ideal spatial adaptation by wavelet shrinkage. Biometrika. 1994;81(3):425–55.
 40.
Anscombe FJ. The transformation of poisson, binomial and negativebinomial data. Biometrika. 1948;35(3/4):246–54.
 41.
Makitalo M, Foi A. Optimal inversion of the Anscombe transformation in lowcount poisson image denoising. IEEE Trans Image Process. 2011;20(1):99–109.
 42.
National Electrical Manufacturers Association. NEMA standard publication NU 42008: performance measurements of small animal positron emission tomographs. Rosslyn, VA: National Electrical Manufacturers Association; 2008.
Authors' contributions
HBL contributed to the derivation and implementation of the algorithm. KW has completed the simulation and experiment together with HBL to verify the performance of proposed method. JT has provided a experimental platform, and made sure the reliablility of experiment data. All authors read and approved the final manuscript.
Acknowledgements
Authors thank Dr. Karen Milada Von Deneen in Xidian University for proofreading this manuscript.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The raw data used for analysis to draw the conclusion has been provided in the Additional file 1.
Ethics approval
The study was approved by the Animal Care and Use Committee of the School of Life Science and Technology, Xidian University.
Funding
This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant No. 2011CB707702, No. 2015CB755500, and the National Natural Science Foundation of China under Grant Nos. 81227901, 81090272, 81527805, 61231004, and 61401462, and the Scientific Research and Equipment Development Project of Chinese Academy of Sciences under Grant No. YZ201359.
Author information
Additional file
12938_2016_221_MOESM1_ESM.rar
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Positron emission tomography (PET)
 Higherorder singular value decomposition (HOSVD)
 Anscombe variancestabilizing transform (AVST)
 Similarity
 Weighted
 Denoising